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Preface 


Volume 1 entitled, "Contamination of Megavoltage X-Ray 
Beams By Electrons And Scattered Photons", was originally 
written as a thesis for the partial fulfilment of a Master 
of Science degree. At the defense of the thesis in 
September of 1982, it was felt by the examining committee 
that Volume 1 should be used towards a Doctor of Philosophy 
degree. 

Section 3, in Volume 1, addressed the problem of the 
transport of charged particles in the build-up region of a 
unit-density phantom. The committee felt that this Section 
provided a good foundation for a more complete study of the 
transport of charged particles generated by photon beams in 
an inhomogeneous medium. A decision was made to undertake 
this investigation*® to “complete ‘the’ thesis. Volume a 
contains the results of this study as well as a method to 
calculate the primary and scattered dose in heterogeneous 
media. 

It was decided to break the thesis into two Volumes 
because Volume 1 0g a self-contained study of photon beam 
contamination. The contents of Volume 2 are not’ referenced 
in Volume 1. For this reason, the Discussion, Conclusions 
and Appendices concerning contamination are left in Volume 
Lee However, in a number of instances, the reader of Volume 
2 is referred to Volume 1, therefore, the pagination in 


Volume 2 is a continuation of Volume lhe 


xXxv1 


peorrerere ton we © omirlo? to edge no. aT -Hoit Gakiasg00 j 


vil@ntgive skw ."s aeelnee bene ta ef ei $ Pr 


et 4 d SIS ong Io eaiteteb ‘edt rr as ; 
tipgoo goldimexns (addy ee i592 eae a1 Bg 
fAweolint 9 sopon€ - wy atiepwet mr Kt ‘ee bi ort 


ih ar ae * ssbbliees a wolfe? oe it ididnge 


& lo isigat di bhi of2 \e8F eon i tueg racneenel io, jnoqaead a 


= 


yolipae aide dade 263 epdedumed eAt | cieltsiato, aan ae 
és Me (here siptqhon simm 3 wt emencesh ‘id a ene a 


eaatsstidd 34 stax ee io mike: axon ane ‘coca 
S: ‘ermaioy oteed?. wet a ms eer : | it 


ot holiten 8° cect Tae Uvcnsatos ‘ee Eyeet) srit aud exes 4 

anisufagpce ton Al +b aired ety Co taint as - 

| a oi | -Siben "| 
€ 7 at [ 

eens {OV ond oiai atieai’ outs —- oo beblven ‘esw f1 


mee sosotte | ‘ae “bore boat ntaened tee pei ! depeitio’ Sah Boo |; 


~ ghoten tne’ nblimusnie adt eosin ahs wt 1 emp ker ob 7 ; 
arate LCF te Pst one 401 amtmadini rors aetusgecth + nu : 


. 


Sempto® te sober, edt sesdnatend 26 aestimien 5 a nt sda 7 
TA seas a suceeilt = ‘ptt fo? os ejashowes ¢ 


Chapter 


5. 


INTRO 


MODEL 


A CON 


ak 
122 


DISCU 
CONVO 


Sed 
8.2 


TABLE OF. CONTENTS FOR VOLUME 2 


DUCTION TO RADIOTHERAPY DOSE COMPUTATION .. 


Sika UEMe Mite On pel UUT POSEY boxe <vefolerere ts lorete: 6 oele eee era's 
BLOC ULONmes Gaia sl.) Dra UWS BEM ise ie terete i. © 0 oe wie 908 
COMCOURMCOMTOC IONS. ‘sierer ord lecie cee wrslet oreisisie. 626-6 
Corrections For Tissue Heterogeneity ...... 


LING DOSE USING THE EGS MONTE CARLO CODE .. 


VOLUTION METHOD -OF CALCULATING DOSE: \.s/<.6:s-0.¢ 


Introduction eeec0e3<e7e0e5ee0nee3es#ee#eeegsc5§c5jeeeeerskeeesseek eee 
Pramamy DOSen oOpread SALT SY Sit cis 610.6°6. 6 tous) ole e ecelone 


7.2.1 Definition of a Primary Dose Spread 
A TL 2) ae edomatalteros svoter siaiter aleve aie erecsrenevetel eheuel eras 
7.2.2 The Generation of Dose Spread Arrays 
Using the MOCA Monte Carlo Code .... 
ite are sR OS ULI iG eer Weber ctalisl store slave! ele: eis) eet oho one olslelece 


Sea teri WOSerOPres Gd rAPTRAV Si 2) oisse.s 0. cisxe evel ee ones 
Convolution Dose Calculation in a 
HOMOPeSNCOU'S GPHAMCOMF Kieren ase %e tele) s)olerelielieteleielcirecevene 
Extension to Heterogeneous Media ......ee0. 
The Spatial Invariance of the Dose Spread 
AIETT AVIS al otiel eo teliei/n/ oh oFehel let evotcha eter ete one) 6] clleiel el olelerelens. sles 
Comparison with Other Methods and 

Porventaad AIMprovVememitsiurs o slenscicle sie) ersielte.s! sa) es\6 
Dose in a Non-Water-Like Heterogeneous 

PIG EOMS voile ators: stole) obs agee! sielal eevee) slates le) 410 leisl ote lense 


SSION AND CONCLUSIONS CONCERNING THE 
LUTION METHOD eeee5«eoee#ers#eeseeos+eoeeeneer5nrleeevw#w##ee#eee ese 


Discussion SOOSOSOGOSOOODUOOOOD OOH COB SOD Oot 
Conclusions eeoeoeoeveevoeveeeeeeeeeeeeeeeeeeeeeeeee 


REFERENCES ecooeoeeoevoeeeweeeeeeeeeeeeeeeeeeeeeeeeeeee 


APPEN 
us 


1 
1 


VITA 


DIGES eee0eee0e#eeeew#enreeeeweweeneenerseeerteeeeeeeteeeeee 


O. Listing of the MOCA Monte Carlo Code .. 
Paes ting. Of; theo program. VolVe's KOK) cstss ss 
2. Comparison of Heterogeneous Dose 
Spread Arrays Calculated by the 
Convolution Method and by the Monte 
CarsVOMMC GUOG: 21016 isis esse. bie ne 0) s\e0le\o.s.6 6/016 61% 


eeeeeere354«eseoecse73e#ee#e#ee#rere#eees+oeteenteeeeneteeteneteeteneteenrenetewesernreeeee#ee#t @ 


OVE 


460 


x 


‘ay 


te 


~ 
- 


yk 


i 
oh) 


a) 


yy ae? 


bene 


2 


ww 


2 sok peo SMEs ahd Aled 6 Nee ae 


OL T ACTED 
o* ve eee eae Oe Ee ; 
ke a ORR aed ae eee ae 


4 , 
sere (es swabs ly eae 


= "2 Lsitsasites.ar oes | 
S200 CBO arKow eos aires 
Lyre Bee SME soa bisa 


> ig pea bs Se ee fer gta 
*-+ ewe ae ee se ee) set 3 ‘bes 


Ls rare aaod as) iat Pe ei 


aa o & 
pn 


La a 


Ld “pee web aks 


ta ‘paige $060 Osi 
.. sbe? ig htad esnem. 


65h fort 
+ hep cei ee ae 


oe he tee? mr 
bésige Srote 
ae © 400 418 Ogee maim 
bow ae 
{6 8 8 Oe eee fee <. oe 
ep oem: agor ibe y 


“a eaten 


Perch. Me ne ae 148 ; 
oe HERI: see hr ea gaa agosalianco) &.8 - | se 


fr a usidnane tials aed WA ome a aanvaie Rae a 
fy ! 

DS 0 s+ 9 Aa art ay A gaDlevdss § 
J 

At : 


Vals 
‘ $4 
oa | 


tae er ete ee © © Be Awe + Ove «a aD 4 8 wed we 46 ee 6 t Pe 


LIST OF TABLES FOR VOLUME 2 


Table Page 


20. Parameters for determining the effective 
attenuation coefficient for 15 MV x-rayS cecccecooes £15) 3) 


xxviii 


¢ FUR We 
GON FO Bade ae 
sIBAT cies | 


ig 


ae 


oe 


i’ 
To | rn 


a 


Figure 


82. 


83. 


84. 


85. 


86. 


Bilis 


88. 


SO. 


90%. 


Cis: 


92. 


93% 


94. 


Jai. 


2G. 


ot. 


LIST OF FIGURES FOR VOLUME 2 


The state of longitudinal equilibrium in a 
ReELOL OS CMeOWUSH MAN TONY coves chojoreue: «peijessi sie4slge oyonepeysideleieiouee/ 6 0. 


Lateral electronic equilibrium in a homogeneous 


phantom eeoeoeoeoeeoeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeoeoeee ee 


Lateral electronic equilibrium in a heterogeneous 


phantom eoeoeeoeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee 


The experimental set-up to measure the 
INHOMOPENE LLVu CORTECLLON. CaACStOLr 2,0.0:0se,0.0 «,04erele eseie.0 6 0 0 0 


The correction factor, for the experimental set-up 
Bel sleal Saat Ae erly ads UT Cau Ok tel oredensdcdensteue: shecenehecsip .supes! ecers.e) 01 8s 


Geometry to take into account contour corrections . 


The definition of the zero-area TAR and the 
secatter-air-ratio eeoeeeeees85§ese#eceoe#sex5#e#25noeees#5ceees858§creeeeeeek ee 8 @ @ 


The geometrical parameters required for the 
Batho or power-law method ceceveccccccccccvcvevevace 


Anew larsi.radssOn. Of .O,WCONNOT: Sorc iGO TEM inccocs s.c.eelec 0s ele 6 


In the equivalent TAR method, the 3-dimensional 
electron density information is "coalesced" 
LIES one cross-sectional slice eeoeoeeoeveeeeeeeeeeeeee e 


The measurement of the perturbation of dose by a 
polystyrene annulus placed in a water tank ......... 


The dose at point P gets a dose contribution from 
an annulus by "revolving" a pencil dose 
QUStri bution about; DOIN Pie .)s sete sco ele co cele ses oosies 


The measured TMR is compared to the EGS 
calculation of the TMR eeeee483+eees83ses#keeg#s#5ceeereeeee#eeeeegeegeree 


A comparison between pencil beam isodose curves 
in homogeneous and heterogeneous phantoms ......... 


A comparison of measured inhomogeneity correction 
factors and ones calculated using the EGS Monte 
Carlo code 1 CLO ene 6. 6. @ 6: © 6: 6 6 0 O ©: 0.6 © OC 0 10 (6 6 OO (0.6 0 O 0 O10 OC 6 010 6 6 (0) 6 


A comparison between the calculated profiles of 
KERMA and dose for a heterogeneous and a 
HOMOgeneouS PNANtOM weeeccceeeeeveeecssesesssseveces 


2G, b,€ 


Page 


263 


266 


268 


270 


Zt 


274 


278 


283 


284 


288 


292 


298 


299 


Oy 


304 


305 


A ne biodiPbaghe te Ty Set. 
San : uis acne As wont whee agin te ea it rey p 
=voehaiuhOe = Mb Saanhine vince somite Tiles 


EA Ue he’ Geet, ga eee Blais a ie eee eee 


avogneyotes edt a ak mada Okage. pave Jatt va 


Tix is. 6. ara ore . Seee eee eee ee eee 


alt D3 eyepnem oF 
KN 09 9 oe me) Pale erg me XOI987 sos 


oli tee Certemiweqae adh me cxaoat sash eee Fe 


Tare rr tee? wll) ee PE eee 
MTS 1. Beka sane andes TAREE DA: on 8 ‘ 
4 ei 
ey poe BAT oer an ote ais 


my FA Ge 8 Ne SRNR 68 Wm de «Nitra eet en ol 


ett wa hart 
ee tinea ies tes Secale SO , 
‘ Tatu 


a 


Wey | 
“Et L + a . a if 4 ya A a8 
Paotreae path ~=) ath) PB it in Bt git 


Pea) 
i) 
> 


a) 


gc ‘eel ny 
“52 Ck stone ak gine ahaha gal ate age 
» | 


‘ 
mond nei >ricievaldidin 4 £ cae ¢ pares a8 Bech. ctt ee 
ged fhwned # “gdiviavas" yd eatuirts ow | 
Res <i te TRE ence ae so tHetis woisydhes j=eis 


sine aak att of teen” et SUT soineash ott Lee 


eee Pre re any to aoFtving?h ss 5 
eattd a maod [ hae reswisd niet ‘~-eemup 2 ae 4 
SOE eens ge Beast 2Hesenege1eisd bos PuQeNS a 65d : 


agisegages Ut ears betnpeon: 4 aGad4 Agmas A oe 
/ POM BOG ady grien beteivoles sano bre erotgn} 


_ ene ieee Hee EAE ORNS Dee eww Ney COO S be eee we Oe a 4 abon olga. 


: to eof itox, tintwiioixs SAT DeasWiad toe +RgmD A TR 


8 bes @¥pMHSZO TSE & jo sob brs AMBRE 


Ge Se Me, Ft Oe Oe E Mie eF 0 ieewlse © meth abe Shoo is BOM $e 


Figure 


98. 


Oo. 


100. 


LG. 


102. 


103. 


104. 


105. 


106. 


LOZ. 


108. 


109°. 


List (OF SE LGURES, FOR VOLUME} 2" (CONT *D) 


A comparison between dose and KERMA correction 
factors for the 15 MV spectrum and its spectral 
components eeee3seoese3<esc§wtsweeeseee37#wexe34nexeee3nre3e+ee3evseeneeeeeteeeeet eee © @ 6 


Primary photons interact in the interaction 
voxel and the charged particles set in motion 
are. Collowed) Chrous tthe? Pham TOM %.. wwetexetetetetereteteteters 


Flow “Ghair't’ jfor ‘photon ttranspor t? part of? MOCAR..5.. 


Flow chart for the part of MOCA dealing with the 
type of photon interaction eeoeoeteee#eoeee#eege#g§5e#een85neeeese#eeee¢ 


Flow chart for the charged particle transport 
part of MOCA eeee7e5eee#scs3ueoeeeeeeee#e#ee#s#eeseeee#eor8#eeeeteeeeeee @ 


The Compton differential cross-section as 
PUTCO) Otek TIVO B1iCePCTION Bye Pesete vito cverebe ele store ere: ole siete. = 


The choice of the kinetic energy for a Compton 
recoil electron eee0e545+eeseeeenseeeeeseetseenrmseeefeenwrseeeeeeeee 


The choice of the kinetic energy for a pair 
production charged particle set in motion ....... 


A comparison of components of the percent dose 
predicted by MOCA and EGS for Cobalt-60 photons . 


A comparison of components of the percent dose 
predicted by MOCA and EGS for 5 MeV photons ..... 


Primary dose spread arrays for 15 MV photons 
for gravimetric densities of 1.0, 0.8, and 0.6 
Whence EVOKE! Gime SO Mats CCI ieee) terete, oe + 0100 6.61 6 


Primary dose spread arrays for 15 MV photons 
for gravimetric densities of 0.4 and 0.2 when 
thesvox ell di mens 1 ONeeiSttL Cm! 27.1.8. a tete oie te ahetete oo Fete Wile" % 


Primary dose spread arrays for 6 MV and Co-60 


photons 66. © SO 6 6 6.6 0) .6 0) 0 6) Oe 01:0) C1916) O.. 8) 0) OO) 'O: (01, 0) OO: 09 0). 0) OD Oe OO eRe 


Primary dose spread array for 15 MV photons 
for a gravimetric density of 0.2 when the voxel 
dimension is 5em eeecoecoevonoeveeveeeeeeeeeeeeeeeeee eee ee 


Truncated first scatter (TFS) dose spread array 
for 1S MV photons Oo. 6 6. 6, 016 © 6 2 0 0 O 0) (6. O'S CO 0 (0 610° 6 0 62) 0) 0 6107 60 


Page 


307 


SZ 


314 


315 


316 


pee, 


320 


322 


323 


324 


325 


328 


35:0 


334 


cna Ye rE a 
- ; ~ ’ at 
(a endo) Ss SOR HOR aaa UOry 


“o>” Dex 


it 
Pia gon AMR AM Ane 5EOb) doowrall ted alnanc A 
; 2 vliogae VN i, SAT, tolmaras of 
Ob, : a SF Bla aes ee ety ee Saito. 
do tse Sin fond Pieces ant {gine veamiat | 
; pt ni 19% Ztent Das Bouentoy sith pine” Tene! 
es » wie KG owl Ay ow Oe Oi 4 erie Tigo" ait AG tt awe Uh eae” A 


to-J1 BG). hOWSS 8ts Metenq 40% 35 wipe wots , } 


usw gtileate FOG Vio” lta? ats. toe sons. eer Od 
sa’ dente ae el a ly tee ee oe? Pe hid ait P 
. 7 eels < Hie Pi ite 


sSeziaiy + $2200 bega8 (7s of ji™od eae wo! he 


G sv eleis spe ee eS ok ee Bagsachy ig Rast E 
& Nejpease-eaet eo tees HI sath bn ae Tenn eh BOS 
Bi < Cusy eae mi Vike eee in sieiw to nota, ine, 

Pas liven ‘ 

d } ) 7 : : iy 7 Shes A : pts 
connie & 743 Pree io solq S) eRe _ 

oa > oa w bets 6 sles ee i. « Bi ea erie ‘4 natal Ware a: 


y Ava 


vl pom tape Det bch aers 4 aid | ‘aot be ait Ah 
GOOF vee bee aol OE Ram ae ta Biciass piles iy) 


+ ae ies 3 
Kae ins o4 90. add Ie Sipegoamen. ae oe lta A BBE % 
SOS sees Shots Ss ealeriiclin Cue ree ALD) Spade : a 
ay 1 we "id ne Wi. 
x Ohi ceeds 2A IOMO. aii tating A i a 
oN a et aig Veale Fa s ‘Tae XA bs Ioibony 
: lee 7 - ae 


1D as 

; | ear Agi 
ia 7 cc) 0) Bak <e 7 7 
buaaiede pal. a . 
: is cz 


7 
7 


a 


. 
ney e.* 
Tie 


ete igt Am 
uke oy oat fone 


Figure 


i238. 


120. 


Wy? 


122. 


122 


124, 


125. 


E26, 


127s 


LIST OF FIGURES FOR VOLUME 2 (CONT'D) 


Residual first and multiple scatter (RFMS) dose 
SspReeadtarrayi sorado MVe phot ONS1 bz fey cress oie Sue:6 os eo slelelee 


Illustration of the dose contribution from the 
PUTeRAe TERA POlR Sp Of 4d Ma CWry ce 0 fo tote so s08 tata oo 050 joe) ose ov 80 e 


The variation of the effective attenuation 
Coceii cientaassaad unetiona Of + dept Ds te oteteus% oitiate oso e 0.60 


The variation of the measured effective attenuation 
coefficients as a function of nominal beam energy .. 


Illustration of the dose contribution from the 
the-aoce (deposi ELON point) Of ViVCW, cies. «1 evctereis seis. syererels 


Measured and calculated TMR's for a 15 MV beam 
as-actunction) of, depth atone the central axis” «dc 


Measured and calculated percent depth-dose data 
for a 6 MV beam as a function of depth along 
the central axis eeee3eeo3s+oeteer#enete#3eesectieeneers35veees#eoee3srn5reeseeeeeee @ 


Measured and calculated dose profiles at dmax 
leone eh 1S) MV beam eoeoee#eeee#8ceeeeeeegeee#+5es#vsnsoeeseseeteeeeteeeee 


Dose profiles at dmax in homogeneous phantoms 
Wit havari Ouse. GENSAt 1 OSes cys ogetecode wrsjeteile crete eiegeyele e066 66.0 6 


The relative primary fluence profile for 15 MV 
beam wedges and the fluence profile used to 
obtain the calculated isodose curve in Figure 123 .. 


The calculated isodose curve for the fluence 
profile shown in Figure 122 and the measured 

isodose curve for a 60 degree wedge for a 

LSM CDCl ivers aleve ee ee © o-ele lots) tualelete: shale oe elels' ele cseretetslete sisrcre 


The measured dose profile at a depth of 5 cm 


when a shield is placed in the 15 MV beam .......... 


The bar shield represented by three fields ......... 


Components of the percent depth dose in a 
homogeneous water phantom predicted by MOCA 


fora 6 MV beam 67610 107 OLE. <6:.8 OLS (0. TELOTO?O © [CO (OF OTC OF OV OFS (OTe 7670 107 e 8 Fe ©. es © 


Components of the percent depth dose in a 
heterogeneous phantom predicted by MOCA 
for a 6 MV beam Sle) 661 Os O10) OO 10 O07 O10). CO O10): 6) 01 O7.0) COTO 8) CFO NOTA) Oh C7 Oe. One. 


XO 


328 


aed 


342 


346 


349 


350 


SO 


352 


354 


355 


SION 


359 


SOW 


362 


ate 


SEs 


4 # © 4.6 


saad Vile Stk ee af rs 
etee nah ilaia atte oF sds 
<7 


ah oenbows ab” tng 


, 
o «Ga 
s ¥ ' 

ve 
e « 


4 


m pees eatp on 


Amt: d reBtigy 2 1 at 


ee a 
9 


(PIB SV ae ai 
Vea msoo nea a 


7 ier oy, 
ro pia. he 44, yee 
4 49,2 4 : | + ¥v AY, 


noi ay Et gqeb 6 nek: 


eid ene lee hd Ae 


+4 Sah a bonnet 


Onar 


oma 


Z ie 


Pty 


Vi ais ein ee : 


ol bus 


uae iF ita 


i « 


rir pr Eine 


ae 


es pina 0 


ia ae aera et aniline Ri eitie aad eff Oe 


— As afrenoginod a a 
Aegomed 


i Roe ‘ “a Be Pe ee ee p J 8 % 1x02 


g of ieee ah airies oh? Jo etaehaqm? . TS! 
RECHE haTaibS ey Mosirade AkgaisyoRss op 


es wig ea eT mle nee eee tee es ore RE ve 8 eo t02 


LIST OF FIGURES: FOR VOLUME 2. (CONT "D) 
Figure Page 


128. Determination of the average density between 
the interaction and dose deposition Voxels: ...c<«ce e's: GO4 


129. Ray tracing is performed by sampling the density 
between the interaction and dose deposition 
voxels eeee<eet3e+oeetee3e#eess3coeeeeeee3nreeee3uaoeoeseeeeeeeeeeeteeeeee 366 


130. Schematic representation of the slab phantoms 
tested to verify the approximations used in 
determining the dose in heterogeneous phantoms ..... 370 


131. The experimental and measured TMR correction 
factor. LOr.a owe OCI. he OCM) DCAM \'eer6 o14 16. oe eset slereisvele-cl oho 


132. The experimental and measured TMR correction 
faecron ior a: LoouMV Ocmeaxm LOCMAb CaM 6 eisrs 0 of¢: s/o) euokerelerelesy oon 


133. The calculated homogeneous and heterogeneous 
TMR profile icone Sy Ike MV eeoeoee0e8es30eo0ese30e0ee3e3weo+oeee3#w#e3nreeeeeeeee AO 


134. A comparison of separating the task of determining 
the wedge dose distribution and calculating an 
inhomogeneity correction factor and performing the 
SpenCila cone ll wOn Ga-Sice prs slersveleieis ia ave ce) sso stensiels oe ereterencrss <O! OO 


135. The dose spread array should be tilted to simulate 
the dose deposition from a divergent primary 
pencil beam interacting (ate 255. se. 2 ss seine siete eee «Oe 


136. A comparison of the convolution method with 
existing dose calculating methodS .....eececvseseees 386 


137. A flow chart of the host computer calculations 
required if an array processor is used to perform 
part of the convolution calculations .«..e.e+creesees 390 


138. A flow chart of the array processor calculations ... 391 


139. The part of an average density array calculated 
for one beam, that intersects with the 
average density array of another beam, may be 
reused by the second beam ..cececcrcscescsscsesssesss JIA 


140. The flow chart of the determination of the dose 
spread array and primary attenuation coefficient 
based on the transmission spectrum of a linear 
BG COTS t OT ies ete oinl cr aroraies ss lols a' 9 enous Stele oie tevele evenetalslei sien sveta Full 


ASO s 


ds 0 


Po fees 0! We 


/ 
" 7 ay 
‘. . ~ oo é ' ‘ i 7 
: WSyted Yt beneh eevee avs Stig Ya oc ad 
‘ , 4 } OTE SORR Bee: FOLRS 
hy Sth. at. Gites ue pan bese at ‘nei om 
(ob) sae ae no Pode sai (bie isew fe ii 
«ws ¢@ © 8 «ee e.« * s«4 . te 
SAGs cule yent Fo ‘Hobe 1: SBaTe aor \ 
ni poet hod) aa Ferree eatd Yat rew | ot A i287" 
ae os Seen Moat Lae “Ags “eQrat aris TA sro, ait antes tm canted 
Wi2o59ias MNaa s See bem bins LRAaea) Leneneeies 
PE vwthiacdkene eee” ae oe ae. ‘Mes a At ed DRI 
= y rs) , 
Ht Sori af potue SSM, fshe PRIRon I ey oT . 
iy haa ple, Pee BS ST RoE tae BOO IM Glos (On uF roy 
sOsSKS7HTSI HH ba <rleetaron ToT AN bor sitet an aT. 
pe 7a 9 # © © 2 eo > @ BP me ** eee 5, VM Gi nia: s! i e71 Pe 
6 : ‘ ay 7 4 ade Pe La 
gh Ledwis29n 2O A884 Sait, seties ies tcvabe twp mode & . Se 
OA BHALPE oe Pes, be Hob puss 1252 tx toile 94098 eat, " 
od, Daleigad + te.) Diem ek “fhe, brgytt: peng s iS Beer i) § 
es ** ees a Wine bs np ee 3 my his ye a 4) 
, ae | > se 
re lumicr od feet Sa: nino mcd an ae ere 
CaM Sy rhe V4 ve Mord ndis teagoh) ezohe 44h, 
Se “see 26 04 = hye — ye ni 47? < iy age ts Saital nt Pboae 
: TE | ie } on Le 
| afew) Oe Obst gat | bee Eis aia deat 
€ Aeraewp se ae a aw ee « ED ee ae seats ie he Ds abide 
hae atid tht Lariat 12% somo 
| fattaSa, OT Mee 


BRR ns tee ao ahotin ore oe : 
2 iy he . cyHnt . : 
DSS eve guekiet on! Fas qoOzeIo07 ama ont tots 
roe | Wiel wa ee Vig 
pstaluatay ANE, pehen, gs gto a ag 1) 
| az at? ites “2Joazis eis ie 

30 Xx ser , Ber ty pulToas y yi lene 


VOLUME 2 


A CONVOLUTION METHOD OF CALCULATING RADIATION 


DOSE FOR MEGAVOLTAGE X-RAY BEAMS 
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5. Introduction To Radiotherapy Dose Computation 


Wetenoed!'<l13-DiICr mror-7 CAL. RTs 


(3-Dimensional Computer Tomography for Cancer Radiotherapy ) 


L.E. Reinstein (62) 
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5.1 Statement of Purpose 


The goal of radiotherapy planning is to obtain a high 


therapeutic ratio, defined as: 


Tumor Control 


Therapeutic Ratio = yooay Tissue Complication 


GS21es) 


This can be achieved by giving a sufficiently high and 
uniform radiation dose to the target volume and a lower dose 
to the surrounding normal tissue. This implies localizing 
the tissues to be irradiated accurately and optimizing the 
dose distribution. 

X-ray images have traditionally been the main technique 
for the localization of tumor tissue. They are obtained by 
passing an x-ray beam through a patient and detecting the 
transmitted x-rays in a 2-dimensional plane perpendicular to 
the beam. The x-ray image represents the transmission 
through all the tissues where the x-rays are passing. AS a 
result, the tissue volume through which the diagnostic x-ray 
beam is passing is "compressed onto a 2-dimensional view. 

A radiotherapy "simulator" is a specialized 
radiographic device which produces radiological images and 
also has the same degrees of freedom of movement as 
radiotherapy units. It? isyprimarily, used, to; localize the 
tumor site with respect to the radiotherapy beam. Even when 
the tumor volume is accurately localized from several such 
views, the 3-dimensional information is rarely used directly 


in treatment planning. Instead a CON TOUL Olea 
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2-dimensional cross-sectional slice through the region of 
interest is usually obtained. This is usually done using a 
mechanical,or,optical),distance...indicator, attached,,to. the 
Simulator gantry which can rotate about the patient. 
Section 5.3 investigates existing methods of correcting dose 
distributions for patient-specific surface contours. 

Internal outines of heterogeneous regions can be 
included by noting the maximum extent of the region in 
lateral and anterior transverse views. The cross-sectional 
view of the heterogeneity is then interpolated by using a 
standard anatomy atlas. The density of such internal 
structures is noted when they are greatly different from 
that of water (e.g. lung or _ bone). The correction for 
tissue heterogeneity correctly takes into account only the 
primary fluence, estimates are made of the effects on _ the 
scattered radiation dose. Because of the uncertainties in 
their location, there was often no correction at all for 
internal heterogeneous tissue contours. Therefore, until 
the recent advances in imaging, there has been little 
incentive to improve the accuracy of the calculations of 
dose deposited in heterogeneous tissue. 

Imaging techniques such as computed tomography (CT 
scanning) (55,56), nuclear medicine tomographic imaging and 
nuclear magnetic resonance (NMR) imaging have improved _ the 
localization of tumor volumes. This alone is expected to 
improved the success of radiotherapy. This information can 
also be used to simplify three aspects of treatment 


planning. 
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1) The patient-specific anatomy can replace estimates based 


on cross-sectional anatomy atlases. 


2) It can eliminate the need for measuring and modifying the 


external contour. 


3) It can eliminate the need to assume the density of 
internal tissues. X-ray CT provides "pixel-by-pixel" (a 
pixel is a "picture element") information on tissue density. 
The problem of extracting the density information from CT 
scans has been studied by a great number of authors (57-61) 

Tomography also provides information in 3-dimensions. 
Therefore there has been an incentive to use this 
3-dimensional information more fully in treatment planning 
(62). In the discussion of dose calculation methods 
(Section 5.4), this aspect will be emphasized. 

Another important advance in radiotherapy, which can 
improve the therapeutic ratio, is the increased availability 
of reliable high energy linear accelerators. A greater 
amount of dose from these sources is deposited at greater 
depths in the patient. Therefore, deep-seated tumors will 
receive a greater dose. At the same time, the build-up 
region extends to a greater depth reducing the dose near the 
patient surface. At higher energies, the range of charged 
particles set in motion can be several centimeters. The 
problem of charged particle transport will be introduced in 
Section 5.2 and investigated using the Monte Carlo method in 
Section 6. A dose computation procedure to take this into 


account will be proposed in Section 7. The effect of 
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charged particle transport is greater for 15 MV x-rays than 
for 6 MV x-rays, and therefore, this investigation has 


concentrated on 15 MV x-rays. 
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5.2 Electronic Equilibrium 


The energy transport in an x-ray beam is a_ two-stage 
process. Photons must first interact to set in motion fast 
charged particles. The charged particles then interact 
continuously with matter and deposit their energy as dose. 
This process was investigated in Volume 1 to explain the 
increase in central axis dose in the build-up region. It 
was assumed that the charged particle transport was 
predominantiy; in the forward direction. iWhe justification 
of the assumption was as follows: If the distance from the 
point of measurement to the field boundary is larger than 
the lateral range of charged particles set-in-motion, then 
for every charged particle leaving the region of the central 
axis there would be a charged particle arriving from lateral 
regions. When this condition is true, lateral electronic 
equilibrium exists (63). In Volume 1, the forward transport 
of charged particle fluence from the point of production was 
approximated by an exponential term. The constant in the 


+ 
exponential, Fe ee describing forward transport, was many 


times larger then the coefficient, a » describing the 
primary photon attenuation. 

Equation iwo<.Secan ibe usedato ‘estimate the “state Jot 
longitudinal electronic equilibrium in a heterogeneous 
region (illustrated in Figure 82). It will be assumed that 
all regions of the phantom are composed of material with the 


same atomic number. Therefore, there is only a difference 


in gravimetric density. The amount of primary energy 
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The state of longitudinal equilibrium in a 


heterogeneous phantom. 


Figure 82. 
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fluence at a depth, z, inside the heterogeneous phantom, is 


the sum of two contributions, one from a unit density region 
I 
(Region I) of thickness z-B, ¥fe+(z)] oe and one from a 


region (Region II) of lower gravimetric density (compared to 


Te 
water) , ©, and thickness B, ¥fet(z)], 


u (Ww) masta “Hot(Z-B) —-u,# (6B) 
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The term, e , Sine Equation ©9.2.t describes the 
attenuation in the low-density region of the charged 
particle fluence generated in the unit density region. The 


attenuation of the primary photon beam by the unit density 


=U! 3 
region is accounted for by the term, S y(z-B) in Equation 
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ine correction~factor CF' is equalt to the ratio of the 
charged particle fluence in a heterogeneous phantom to the 


charged paticle fluence in a homogeneous phantom: 


I+I] 


+ 
Vie vz] u.0l-9) B 
CF (z) = BETA ee SLD =e iv 
¥le (2 Ion 
Seo) 
Equation 5.2.5 has no dependence on ren ° This means 


that if longitudinal equilibrium has been established, then 
longitudinal electronic equilibrium will not be perturbed in 
the low-density material. Equation 5.2.5 predicts that the 
inhomogeneity correction of the primary dose should always 
be greater than unity. 

Perturbations in lateral equilibrium are not so easily 
described quantitatively. Analytic. solutions of equations 
describing the transport of charged particles have yet to be 
achieved for generalized boundary conditions. Only Monte 
Carlo calculations provide a method to accurately quantify 
the role of lateral equilibrium. However, a qualitative 
description of the establishment of lateral equilibrium in a 
homogeneous’ phantom and its loss in a heterogeneous phantom 
is instructive. 

Figure 83a) illustrates a homogeneous phantom 


irradiated by a non-divergent photon beam with a field size 
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Lateral equilibrium does not exist at the 
central axis when the field size is less 
than twice the lateral range of electrons. 
Lateral equilibrium never exists near the 
field boundary. ; 
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smaller than the lateral charged particle range. The 
charged particles are assumed to be set in motion at a fixed 
angle with respect to the incident photon beam. For the 
sake of demonstration, it is assumed here that they do not 
deviate from their initial direction. The charged particles 
originating near the beam boundary deposit energy at the 
central axis and travel beyond (the dashed lines) to the end 
fT neLr range. If the field size is increased (Figure 
83b)) charged particles originating near the new boundary 
have enough range to reach the central axis increasing the 
dose deposited there. Laterat, ‘electronic equilibrium is 
established when a further increase in field size will not 
result in a further increase in dose because charged 
particles set in motion near. the field boundary (Figure 
83¢))-do not have suéficient—gangesstoewrreach the ‘central 
axis. 

Figure 84a) illustrates the loss of lateral electronic 
equilibrium inside a heterogeneous region with a density of 
one half of the density in Figure 83. Since the range of 
particles in the low-density region will be double that in 
the eunitrevdensity “region, “albiteldmescize. sufficient to 
establish equilibrium in a. unit cdensity medium must be 
doubled in order to obtain equilibrium in the low-density 
region. 

Even though this analysis of electronic equilibrium is 
heuristic it would be expected that the establishment and 
loss of lateral electronic equilibrium in heterogeneous 


Phantoms may. be “a-imore important effect than that of 
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A field size large enough to establish 
lateral equilibrium at the central 
axis may not be large enough in a low- 
density region. b) If the hetero- 
geneous region has a density ot one— 
half, the field size must be doubled 
to re-establish lateral equilibrium. 
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longitudinal equilibrium. 

A series of experiments was used to verify this 
hypothesis. Figure 85 illustrates the experimental set-up. 
An ion chamber (Capintec PR-O6C) without a build-up cap was 
pracea at the) central axis a constant 100 cm from the 
source. The thickness of the cork region was also. kept 
constant at T2.9)-em- = 0.2 cm. The total thickness of 
overlaying material was 15.6 cm+0.4 cm. The position of 
the cork was’ varied with respect to the ion chamber. The 
quantity, A, is the distance from the first interface to the 


probe. The correction factor, CF, was defined to be: 


alae ee Heterogeneous Phantom TMR(z ,A,W) 
Homogeneous Phantom TMR(z,W) 
(Secs 0) 

Figure 86. iVustratest the: .conrectiongwfactor. as a 
function of A for various fieldisizes. At small field sizes 
inside cork, the correction factor is less than unity. This 
means that the dose in the low-density heterogeneous phantom 
is less than the dose in the homogeneous phantom even though 
the photon fluence in the low-density heterogeneous phantom 
is greater than in the homogeneous phantom. A small 
increase in the field size at field sizes less than 
10cm x 10 cm (for example from 5 cm x 5 cm to 6 cm x6 cm) 
results in a relatively large increase in the correction 
factor... At a frela size of 10cm x“ l0.cm~> therev= ts  Tittite 
further increase in the correction factor with field size. 
Both these effects are expected if lateral equilibrium did 


not exist for the smaller fields. 
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Figure 89d. 


The experimental set-up to measure the 
inhomogeneity correction factor. The 
source-to-probe distance, the thickness 
of cork, and the depth of the probe from 
the surface was fixed. The position of 
the cork with respect to the probe is 
varied. 
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Figure 86. The correction factor, for the experi- 
mental set-up illustrated in Figure 85, 


for various field sizes. 
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Similar measurements have been performed by Young and 
Kornelsen for 10 MV x-rays (64). They also report a loss in 
lateral equilibrium in lung-equivalent materials when the 


frel'd size i's’ small. 
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5.0 Contour Corrections 


Measurements of absorbed dose from a megavoltage x-ray 
beam are usually obtained in a water tank or phantom. The 
central axis of the beam is perpendicular to the water 
surface and the dimension of the tank is larger than any 
beam width to be measured. Figure 7S$7a) illustrates ‘the 
geometrical arrangement for the measurement of dose at point 
Poy ine surface ‘of a patient is rarely flat, and thus ‘the 
central axis, where the beam enters, is’ generally not 
perpendicular to the surface. The measured dose must _ be 
corrected to take patient surface contours into account. 
Figure 87b) shows the arrangement for the calculation of 
dose at point P which nows lies under a curved contour. A 
"ray" from the radiation source through P represents a 
pencil beam. 

The "effective attenuation coefficient method" takes 
the “surplus” or “deficit “of ebkissue Pinto VWaccount by uthe 
subtraction or addition of dose) (63). ethe true depth,” Z-; 
of a pencil beam through tissue”> with -a ‘contour is 
determined. The difference between the depth in the 
all=-water situation from a reference surface, Z, and the 
true distance, Zz", is multiplied by an empirically 
determined coefficient, | Oras that depends on the photon 
energy to obtain the contour corrected dose: 

|e Ara) 
DOSE(z~*) = DOSE(z) e eft =3DOS5E (2) [es koge (2-2 0] 
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Water Tank 


Patient Contour 


Figure 87. 4) Measurement of the dose in a water 


tank. b) The dose must be modified 


to take into account the patient sur- 
face contour. 


- ¢ in a ena > oy 
mn Tah PON P 
. : a 7 ve ¥ 
7 i? Shs 
a pe t 
as: a ne 
-¥ he 
: he es 


s 
— 


# 
Seo oem aie 08 re 


ie | 


— 
‘ 

# i. 

al, 4 


tt 


o 


THi#W Pf 2 DBOb 609 26% IAsineee d 42 Le ica 
‘ 5 a x _ y * 4. 
LOFT hem ad Fenn eset mt ue = ial call aa: 


—Wie Troi ipa 


en, any fnlas SH, Os Roe | - _ 13 a 


‘am ue a) ~ee fee 


A 
3 4 _ a 
| j a "sie ae 
> A <2 “ = 
a aa = a ' i ie - _" ‘ sua — 
a id o : 7° - : a - Pan 
¥- oa ie = i @ ° na ¢ - iV 
7 ies a > y 7 F Mar Y 


Ze 


The coefficient, Kee for Co-60 beams is smaller than the 
linear attenuation coefficient so the effective attenuation 
coefficient method must take an unspecified amount of 
forward directed scatter as well as primary attenuation into 
account. The dose must be corrected at all depths for’ each 
pencil beam. 

The "isodose shift method" (63), shifts the depth of 
all isodose lines by an amount directly proportional to the 
difference between the true and Sétetence distance, 2z'-z. 
If the true thickness of tissue is less than the all-water 
thickness, the isodose lines are shifted towards the surface 
and in the opposite situation, the lines are shifted away 
from the surface. The method is computationally simpler 
than the effective attenuation coefficient method because 
the shift does not have to be applied at each depth. The 
constant of proportionality depends mainly on_ the beam 
quality and like the effective attenuation coefficient in 
the previous method it must be determined empirically. The 
isodose shift method is inaccurate very near the surface of 
the phantom. 

The main disadvantage of the effective attenuation 
coefficient method and the isodose shift method is the 
determination of their empirically derived coefficients. In 
general, the coefficients depend not only on the energy of 
the beam but also to some extent on the beam size, depth of 
the” calculation point, and the source-to-surfacedistance 
(SSD). 


Two other methods, the "ratio of TAR" and "effective 
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SSD" methods, also employ the true length through tissue, 


z', but explicitly take into account the beam size, and 
therefore, account indirectly for some scattered radiation. 
The ra two of TAR method is based on measured 


tissue-air-ratio (TAR) deta The tissue-air-ratio under a 
eontour) diisyh-tounds” by fanterpolating: ~ in’) at‘ look-up’ “table, 
eompirled) ) asie ai = functions of* depth’ and’ fiveld’ size, using the 
true depth. Since TAR's are usually measured along. the 
ecentraile axisi,(ethe method is, strictly-speaking, limited to 
correcting the dose _ there. However, in -practice,,® thii's 
limitation is ignored and off-axis corrections are estimated 
using the central axis data. 

The effective SSD method corrects the tabulated percent 
depth-dose data in a manner’ similar to the ratio of TAR 
correction method. The percentage depth-dose must also _ be 
corrected to "remove" the inverse square attenuation. The 
corrected percentage depth-dose P(z',W,SSD') is: 


2 


Pt. WeSSD-) = pe’ .w,ss0y( $$0* 2, ) 


leopersras) 


* TAR measurements are obtained in a similar manner as 
tissue-maximum-ratio (TMR) data (see Figure 11), except the 
normalization value is: ~obtain with the probe in aor 
surrounded by a build-up cap of sufficient thickness to 
establish electronic equilibrium. Like TMR measurements the 


measurements for TAR's are obtained at a fixed distance from 


the source. 
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This correction can be made to off-axis positions as well as 
positions on the central axis. 

Since TAR's and percentage depth-dose data contain 
scatter as well as primary dose, these corrections may 
correct somewhat for perturbation in the scatter dose. The 
major drawback of these methods is a failure to consider 
explicitly the influence of the scatter contribution to the 
contour correction. 

The "differential scatter-air-ratio (dSAR)" method 
takese the contour *configuration into ‘account (65). The 
differential scatter-air-ratios are produced from 
scatter-air-ratio (SAR) data which are in turn derived from 


tissue-air-ratio (TAR) data: 


SAR (257) (= FAR(z,1r) = TAR(zZ;0) 
(52555) 


Where, r, is the field radius. TAR-CZ,,OD 8 LS Teal leds tie 
"zZero-area" tissue-air-ratio. It cannot be measured 
directly but is extrapolated from TAR data as the field 
radius is decreased (see Figure 88). The zero-area TAR 
corresponds to the primary dose if the charged particle 
energy is locally deposited (i.e. if electronic equilibrium 
exists). The SAR values are then a measure of the dose 
contribution of scattered radiation in a beam. 

The dSAR. values at ayspecific. depth,, Zz, ares equal “to 


the difference in SAR values per unit change in the field 


size: 
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TAR(z,W) 


SAR(z,W) 


Field Size W 


x Measured 


ah Extrapolated 


Figure 88. The definition of the zero-area TAR, TAR 
(z,0) and the scatter-air-ratio SAR (z,W). 
TAR (z,0) is found by extrapolated TAR 
measurements to a field size of zero. 
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Equation 5.3.4 represents the scatter contribution for a 
pencil beam at radius, r, and angle, @. The total scatter 
dose is found by summing the contributions as a function of 


radius and angular position: 


2nR eee i 
Decatter 2») = Aa eee —eaieipe ee 


(36.325) 


D. is the dose in air at a distance from the source that 
eorresponds to the depth, Zz. Tce is the relative 
primary fluence at the radius, Ti) and angular position, cae 
incident on such pencil beams. 

One of the main difficulties with the dSAR method is 
the extrapolation of TAR's to get the zero-area TAR values. 


The extrapolation takes place where the TAR values are 


changing rapidly in a non-linear fashion making. the 
procedure technically Gifhrewke. There is a more 
fundamental problem. The TAR concept is defined only if 


electronic equilibrium exists, but this is impossible for a 
pencil beam of zero-area. Electronic equilibrium only 
exists’ in a lateral direction if the field size ais larger 
than the lateral range of charzed particles set in motion. 
The TAR's are extrapolated to field sizes where equilibrium 
does not exist. The dSAR method was first introduced for 


use with Co-60 beams such that the pencil beam is 
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sufficiently wide to establish equilibrium. The assumption 
of -erectronic equrlibrium starts to break down for higher 
x-ray energies. 

A second problem is that central axis data is applied 
to off-axis points. This may not be serious in practice, 
but the approximation has not been validated by modelling or 
experiment. 

A final problem, is that the dSAR method does not’ take 
into account the contour where the beam exits the patient. 


Therefore, backscatter is not dealt with properly. 
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5.4 Correction For Tissue Heterogeneity 


Contour corrections are Simply a special case of 
corrections to account for heterogeneous tissue densities. 
The density outside the patient can be considered to be 
approximately zero. The only fundamental difference between 
contour corrections and other low-density heterogeneity 
corrections is that the dose in air outside the patient is 
irrelevant to the treatment and need not be calculated. For 
this reason, the effective attenuation coefficient, isodose 
shift, TAR ratio, and effective SSD methods have been 
applied to tissue inhomogeneity corrections. The true 
distance in the tissue, z',(discussed in Section 5.3) is now 
interpreted as the "radiological" or water-equivalent 
distance in the heterogeneous region. The major fault of 
these methods, when applied to heterogeneity corrections, is 
a failure to take into account the relative position of the 
heterogeneous region with respect to the calculation point. 

The Batho or power-law method takes the depth of the 
heterogeneity, with respect to the depth of the point of 
measurement, into account. The original Batho method 
assumed that the heterogeneous region was a Slab at right 
angles to the beam with an extent larger than the field 
boundary (66). Electronic equilibrium was assumed not to be 
perturbed and that only the dose beyond the heterogeneous 
region was required. The method was extended to include the 
dose inside, as well as distal to, the heterogeneous region 


by Sontag and Cunningham (67). Their formulation for the 
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correction factor, CF was: 


=) 
TAR(z,,W) * > 
GFize.z5 4) \= = 
Lee TAR (2, ,W)? on 
GSe4e) 
24 and Z> are shown in Figure 89 and Py and P5 are the 


relative electron densities (compared to water) of regions 1 
and 2, respectively. 

The Batho method has been subsequently modified by 
Cassell, Hobday and Parker (61) and Lulu and Bjarngaard (68) 
to take into account multiple slab geometries and estimate 
the” correction: factor ‘when the Fateral! -extent (of the 
heterogeneity is less than the beam width. Limited success 
has been achieved by modifying the Batho method to take into 
account electronic disequilibrium (69). 

The equivalent TAR method developed by Sontag and 
Cunningham (70) was inspired by O'Connors theorem. O'Connor 
(71) proposed that the dose at corresponding points in two 
media with different density, but the same atomic number, 
Wille eber thes same everywhere provided ariel distance 
measurements are scaled with the density (see Figure 90). 
The basis for the theorem is the linear dependence of the 
attenuation coefficient on density. 

The premise on which Sontag and Cunningham based their 
method, was that an. equivalent TAR could be found in 


tabulated TAR tables with the depth and field size suitably 


scaled: 
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The geometrical parameters required for the 
Batho or power-law method. 
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SSD’ 


kw 


SSD’= © SSD 


We 


Dose(z/.r’W’,SSD ) = Dose(z,r,W,SSD) 


Figure 90. An illustration of O'Connor's theorem. Al] 
the distance measurements in this figure are 
inversely proportional to density. 
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1 TAR Ww) 
CF(2"W) = aR 
(5245.7) 


The tissue-air-ratio was divided into primary (zero-area 
TAR) and scatter components: 
TAR(z’,W) = TAR(z~,0) + SAR(z~,7) 
(53453) 


(Se4n4) 


The term, r, isia circular field radius that produces the 
same dose at all depths along the central axis ina 
homogeneous phantom in a rectangular field. The term,é , is 
an effective relative electron density, weighted with 
respect to the scatter dose contributions from tissue 
elements surrounding the calculation point. The weighting 


procedure should be carried out over 3-dimensions: 


gk UK (5.4.5) 


The weighting factors, Wak? is) the Gcontributlong aco, the 
scatter dose from points surrounding the calculation point. 
The weighting factors are large in regions proximal to _ the 
calculation point and for points close to the calculation 


point. 


Weighting in 3-dimensions, to arrive at an effective 
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electron density for each point, was felt to be too time 
consuming for computers available at the time. Since: Cr 
information is produced and stored in cross-sectional 
Slices, the 3-dimensional dose was estimated using this 
plane-by-plane format. The electron density information, 
for all the slices neighboring the slice in which tthe dose 
calculation was being performed, was "coalesced" into a 


2-dimensional effective relative electron density array (see 


Figure 91): 
lg CG pee 
oe 
€.- Se T 
ak » We 


(5.4.6) 
This 2-dimensional array is produced before the dose is 
calculated. Le represents the scatter dose weighted 
electron density as if all the scatter dose was coming from 
one neighboring slice located at an effective location. A 
weighted relative electron density is determined for each 


calculation point in the slice: 


(S.A 7) 


The above procedure is repeated for each point in the 
calculation splane. It has been assumed that the weighting 


factor can be separated into space components: 
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(5.4.8) 
The term W; is given by: 
He = SAR(z=10cm, r,) - SAR(z=10cn, r,) 
($24..,9) 
The radius terms r, and r, are related to the equivalent 


ii Z 


circular field and the distance between neighboring slices 
in a non-linear fashion. The effective distance to the 


coalesced scatter slice, Nace is given by (see Figure 91): 


(Sar 0)) 


The calculation is speeded up by the spatial separation 
of the weighting factor because the summation for each 
calculation point is carried out over two dimensions instead 
of three. However the assumption that the non-calculation 
planes can be replaced by one effective plane placed at an 
effective distance away was never validated by Monte Carlo 
modelling or experiment. The procedure could not be proved 
by a comparison of Equation 5.4.5 with 5.4.7 because the 
3-dimensional weighting factors were never obtained. A 
Simplification of the equivalent TAR method has’ been 
proposed by Tatcher and Palti (72). Equation ’+Ss4—53 “is 
maintained except, z', is the radiological length through 


the heterogeneous region and, W, is replaced by: 
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In the equivalent TAR method, the 3-dimensional 
electron density information is "coalesced" 
into one cross-sectional slice. 


Figure 91. 
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WeeNe 
(5.4.11) 
a Ve eees er iak 
oe 
Va * 222 Vigx Fagx 
ijk 
(5.4.12) 


Where £€ and °i5k refers to the relative density compared 


to water. V and are the volumes of the 


W Vijk 
water-equivalent and heterogeneous parts of the _ phantom, 


respectively. f is@ equal -to oo tor cork and equal to@= 


ijk 


Cleo =) "tor teflon. sts tater claimed that f. is not 


ijk 
sensitive to variations in density! It is also claimed that 


fi ik 


is nearly independent of geometrical position although 
ube is only applied when calculating” “the” “dose “in 
heterogeneous phantoms! 

The delta-volume method was originally suggested by 
Cunningham and Beudoin (73) and refined by Larson and Prasad 
(74) and Wong and Henkelman (75). The method, as formulated 
by Wong and Henkelman, divides the dose into primary, first 
scatter augmented by some second Seatter,= and “resiaual 


mudtipre "scatter “dose. Mathematically, the expression is 


given by: 
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x 
DOSE = PRI # “ 
ieee bef Pik 1 :i9k fosigkeliaik, im 
(5.4.13) 
SAR (B:d,5-r) ser 
= m a : digi porte 
Sn 7 TSAR Czy Spied) + BEE (PS 5K) ] 
m ijk 

(5.4.14) 


where is the density of the volume element (voxel) 


Pijk 
which is the origin of the first or augmented second scatter 


* 
GCose, at point 1) Ks SS foeeeie  CONLCriLDULION OL 9 first 


and augmented second scatter dose (second scatter dose 
contained within an angle of 45° centered about the first 
scatter ray) to the calculation point i,j,k determined for a 


homogeneous water phantom. is the change in the 


tosijk 
primary V~actenuatvon “along the » primary ~ ray ~to the point 


he a el aan is the change in the rirst scatter 


attenuation along the path between the point i,j,k and the 
calculation’ point. The term Si (med ) approximates the 
residual multiple scatter dose in a heterogeneous medium. 
AH jk 
a "small void. ‘at the point i,j,k. It’ is’ determined by the 


is the perturbation in the multiple scatter dose from 


difference between the experimentally measured dose 
perturbation (see Figure 92) and the calculated value using 


AS ° 6 is the overall mean density of the heterogeneous 


* 
ie ky k 
phantom. SAR (d,r) is the residual multiple 


scatter-air-ratio for a homogeneous water phantom and is 
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obtained using SAR values and the augmented first scatter 


values. The determination of AH, . Ys worthy of note. et 


jk 
is obtained by displacing water by low-density polystyrene 
foam constructed in conical rings about the central axis 
CUB). The small perturbation in Ydose\ due to the 
displacement of the ring is measured. ihe: sradius of ‘the 
ring, the depth of the measurement point and the position of 
the ring with respect to the measurement point is varied. 
The experimental set-up is shown in Figure 92. The number 
of measurements involved requires automation of the 
procedure. The first and augmented second scatter dose are 
calculated and their contribution is subtracted from the 
total scatter perturbation measurements’ to arrive at the 
multiple scatter dose perturbation values. 

The A -Volume method was the first "true" 3-dimensional 
dose calculation algorithm. It exactly corrects the first 
scatter dose due to the presence of heterogeneities 
(assuming electronic equilibrium of the scatter photons) and 
estimates the multiple scatter dose correction. 

The delta-volume method has not been optimized for 
speed of calculation. AS many operations are carried out to 
determine the dose corrections for first and augmented 
second scatter dose as for the residual multiple scatter 


dose, even though the contribution of the latter is much 


smaller. Calculation time should be appropriated on the 
basis “of "Mthe importance of the contribution of the 
calculation to the overall result. However, the speed 


limitations may not “be critical. The method is being 
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Figure 92. 
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Polystyrene foam conical annuli were placed 


in a water tank. The perturbation in dose 
compared to the all-water saturation was 
measured. 
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implemented using a dedicated microcomputer. A group 
promoting the A-volume method (Mallinckrodt) is using Very 
Large Scale Integration (VLSI) microelectronic technology to 
develope a "chip" that will do the "ray-tracing" portion of 
the A-volume calculation rapidly (77). 

A major criticism of all the above methods is that. the 
inherent assumptions are not identified explicitly and they 
have not been validated directly by measurements or 
Simulation (eg. the Monte Carlo etnodon Consequently, 
their scope of applicability is ill-defined, and the 
clinical situations in which they are limited are difficult 
to identify a priori. 

The major fault of the equivalent TAR and delta-volume 
methods is an inability to handle situations of electronic 
disequilibrium. Most of the dose arriving at a point is due 
to primary dose which will be most perturbed by changes in 
the state of electronic equilibrium. A related problem with 
these methods is that they are based on the separation of 
primary and scatter dose using the zero-area TAR and SAR's 
obtained empirically. The concept becomes more and more 
untenable at higher photon energies as electronic 
equilibrium becomes more important and the extrapolation of 
TAR's to such fields is incompatible with the equilibrium 
requirement. 

The need for considering non-local energy deposition 
has been documented by several authors (78-81). Some photon 
dose parametrization models have accounted for electronic 


"build-up" longitudinally along the central axis, but none 
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have “rigorously treated the lateral “spread of charged 
particles set-in-motion (40,82) (see Section 3 in Volume 1). 
Young and Kornelsen (64) have taken into account the dose 
reduction caused by a lack of charged particle lateral 


equilibrium using a semi-empirical "loss factor". 
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6. Modelling Dose Using The EGS Monte Carlo Code 


The state-of-the-art of the Monte Carlo 


method in 1952: 


"The electrons or photons were followed 
through successive intervals and their 
fate in passing through a given interval 
was decided by spinning a wheel of chance; 
the fate being read from one of a family 
of curves drawn on a cylinder"..."the 
[cylinder] motor was observed to stop at 
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It has been shown that lateral electronic equilibrium 
does not exist near field boundaries and along the central 
axis in low-density heterogeneities for small fields (see 
Figure 86). The EGS (Electron Gamma Shower) Monte Carlo 
code was used to investigate regions of lateral 
disequilibrium directly by simulating the transport of 
photons and charged particles set in motion. 

EGS was originally developed at Stanford University for 
the study of the cascade of charged particles and photons 
produced by high energy (up to 1 GeV) cosmic rays or photons 
(C33: )%. Rogers, at the National Research Council of Canada, 
has corrected and modified the code to make it applicable to 
the energies encountered in medical physics’ and health 
physics problems (84-86). 

EGS takes into account most of the radiation physics 
required for megavoltage energies. The interactions treated 
for photons. .are..pair production, Compton, and the 
photoelectric effect, and for charged particles, multiple 
scattering (Moliere's theory as _ formulated by Bethe), 
"knock-on" collisions (Moller scattering for electrons and 
Bhabha scattering for positrons), continuous’ slowing down, 
bremsstrahlung and positron annhilation. The only present 
flaw is an inaccuracy (up to a factor of two) in calculating 
correct bremsstrahlung spectra (86). The charged particle 
"cut-off" energy, below which the residual energy can _ be 
deposited "on the spot", can be varied continuously upward 
from 10 keV. 


The program was implemented to produce the dose 
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distribution in homogeneous and heteragencous phantoms from 
primary pencil photon beams. The modelled beam was a 
0.50 cm in radius, monoenergetic parallel beam. The phantom 
was composed of 20 cm of either "water" or a combination of 
"water" and "cork" (water-like in chemical composition but 
with a gravimetric density of 0.25) slabs totalling 20 cm in 
overall thickness. 

The pencil beam dose distribution was also used to 
compose the -central axis dose due to a broad photon beam. 
Figure 93 illustrates the procedure schematically for a 
pencil beam directed into the page. The small circle of 


area, A represents the pencil beam. The dashed lines 


pencL AL: 
represent an annulus in which dose has been deposited from 
the pencil beam. DOSE*® ‘CONTRIBUTION (P) is the dose 
contribution at point P from an annulus of area A 

annulus 
(shown by solid lines): 


AANNULUS 


DOSE CONTRIBUTION (P) = [DOSE] 
ADENCIL ANNULUS 


RO". 1) 


The ratio of the area of the annulus to the area of the 
pencil beam gives the relative amount of energy being 
released in the annulus compared to the pencil beam. 

A widely used "rule of thumb" states that the mean 
energy of x-rays produced in an accelerator is approximately 
one-third of its nominal energy. Therefore, the mean energy 
of a 15 MV x-ray beam should be approximately 5 MeV. Figure 
94 illustrates that the calculated tissue-maximum-ratio as a 


function of depth is within 5 % of the measured values for a 
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Figure 93. The dose at point P gets a dose contribu- 
tion from an annulus by "revolving'' a pencil 
dose distribution about point P. )The value 
of the pencil dose distribution 1s (POSE ius 
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Figure 94. The measured TMR is compared to the EGS 
calculation of the TMR. 
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field size Sou tveiens to s.ae/ 10.0 cm eidiameter field. The 
present version of the program accepts only a monoenergetic 
initial photon energy, so that each energy component in the 
spectrum requires a separate run of the program. The 15 MV 
spectrum used to model the dose in the build-up region (see 
Figure 70) had 36 energy components. Therefore, this would 
require 36 runs of the program to produce the dose results 
for the x-ray spectrum. It was felt that if the dose, as a 
function of depth (see Figure 94), could be modelled with 
good precision with only one energy component (5 MeV), then 
a spectrum containing six components would improve the 
agreement. The energy range between O and 15 MeV was, 
therefore, divided into six components and the mean energy 
of each component determined. The 36-component 15 MV 
spectrum (Figure 70) was used as a starting point. The mean 


energy of the six components were chosen with the mean 


energy of the bins of the 36-component' spectrum, Che ee 
using: 
EO an ; Fam 
ro) os » Fnm 
m 
(6.2) 


Where bok is the fluence of photons in each component of the 


36-component spectrum. The photon energies that resulted 


* The EGS program assumed cylindrical geometry so the 
equivalent circular field was used to model rectangular 


fields (54). 
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were 0.18 MeV, 0.67 MeV, 2.57 MeV, 6.61 MeV, 11.12 MeV and 
13.78 MeV. A determination of the weighting of the spectral 
components, to give the best agreement with measured data, 
was made using the least-squares’ method. Figure 94 
illustrates that the photon fluence, weighted in the 
proportions 0:11:11:11:5:0, improved the agreement of the 
calculated and measured dose to within 0.5 %. 

The effect of a low-density water-like slab 
heterogeneity on pencil and broad beam geometries was 
investigated. The cut-off kinetic energy for the charged 
particle was chosen to be 189 keV (total energy of 700 keV). 
Figure 95 illustrates a comparison between pencil isodose 
curves in a homogeneous (on the right) and a heterogeneous 
phantom (on the left). The 50 % isodose line is not present 
inside the heterogeneous portion of the pencil beam. The 
5 % and lesser isodose lines "bulge" outwards. Thas 
indicates that charged particles are "streaming" outward 
from the pencil beam into the surrounding low-density 
region. At first sight, the 0.1 % isodose line does not 
seem to be a very important contribution to the dose. 


However, A increases with distance from the central 


annulus 
axis. Therefore, Equation 6.1 suggests that a shift in 
isodose lines. to greater radius “will produce a greater 
contribution when summed for a broad beam. 

The TMR inhomogeneity correction factor for a_ broad 
beam was determined by summing pencil beam contributions for 


a heterogeneous and homogeneous phantoms. The thickness and 


density of the heterogeneity was fixed at 8 cm and 
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A comparison between pencil beam isodose 
curves in a homogeneous (right) and hetero- 
geneous phantoms (left). The shaded region 
has a density of 0.25 g/cm>. 
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0.25¢/cm>, respectively. The total thickness of overlying 
material was 15.5 cm (actually, the average dose between 15 
and 16 cm was’ determined). Figure 96 illustrates’ the 
calculated correction— factor. Also shown is the measured 
correction factor from Figure 86 obtained with cork as_ the 
low-density heterogeneity (with a density of 0.30 70.02 
g/cm). There is a 2 % discrepancy between the calculated 
and measured data. Most of the discrepancy is probably due 
to the difference between the actual cork density 
(0.30g¢/em”) used in the measurements and the density used in 
the Monte Carlo simulation (0.25¢/cm). When determining 
the correction factor from measurement there was an 
underlying assumption that the mean stopping power of 
charged particles in a disequilibrium situation was equal to 
the mean stopping power in an equilibrium situation. This 
is equivalent to assuming that Equation 1.5.15 can be 
Simplified to a simple ratio of ionization measurements. In 
contrast, the Monte Carlo correction factor implicitly takes 
into account the particle's changing stopping power as it 
loses energy. It is unlikely that this effect would produce 
a discrepancy of 2 % (53). 

Lateral disequilibrium was investigated directly by 
setting the charged particle cut-off equal to the incident 
photon energy. This forced all of the charged particle 
energy to be deposited "on the spot". In this way, the 
KERMA, rather than the dose, is determined. Figure vw 27, 
illustrates a comparison between the profiles of KERMA and 


dose for a homogeneous and a heterogeneous phantom for a 5.0 
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A comparison of measured inhomogeneity cor- 
rection factors and ones calculated using 


the EGS Monte Carlo code. 


Figure 96. 
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A comparison between the calculated profiles 
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of KERMA and dose for a heterogeneous (left) 


and a homogeneous phantom (right). 
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MeV pencil beam. The profiles were taken at a depth of 
10.5 cm. In the homogeneous case, the parameter, A, was 
4.5 cm. It can be seen that in a homogeneous phantom, KERMA 
equals dose at a lateral distance of 1.5 to 2.0 cm from the 
central axis of the pencil beam. However, dose does not 
equal KERMA in a heterogeneous phantom even at a distance of 
5.0 cm from the central axis of the pencil beam. 

Broad beam KERMA distributions were produced from _ the 
pencil beam KERMA distributions. Figure 98 illustrates a 
comparison between the dose and KERMA correction factors for 
the 15 MV spectrum and its spectral components. The dose 
earrecGion factor at first increasesfas a function of: . field 
radius and then decreases. The KERMA correction factor only 
decreases. At an energy of 0.67 MeV the dose and KERMA 
correction factors are in close agreement at field radii 
greater than 1 cm. At higher energies, the discrepancy 
between dose and KERMA correction facors becomes greater. 
Therefore, at lower energies KERMA is a good approximation 
to dose, but at higher energies, the KERMA equality with 
dose breaks down. 

The measured correction factor is also shown in Figure 
9387 For higher energy beams the dose increases as a 
function of field size, whereas the KERMA always decreases 
as a function of field size. The increasing trend of the 
masured correction factor as a function of depth agrees best 
with the calculated correction factor using non-local energy 


deposition (the dose curves, rather than the KERMA curves). 
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A comparison between dose and KERMA cor- 
rection factors for the 15 MV spectrum 
and its spectral components. 
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7. A Convolution Method Of Calculating Dose 


Knowledge finds differences, but 


understanding seeks similarities. 
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mn ointrodiue tion 


Section 3 in Volume 1 indicated that simulating the 
longitudinal transport of charged particles set-in-motion 
enabled the dose in the build-up region to be _ predicted. 
Section 6 in this Volume illustrated that the lateral 
transport of charged particles must be included to. predict 
lateral disequilibrium phenomena. Therefore, a dose 
deposition model for megavoltage x-ray beams is required 
which explicitly takes into account charged particle 
transport. 

Ine thisiy secttomata method trof tanad ysis iwhich is 
physically sound and internally consistent will be 
discussed, which allows calculation of 3-dimensional dose 
distributions in homogeneous or heterogeneous’ phantoms 
frradiatede with ‘stilted sesofi t “aigh energy photons when 
electronic equilibrium is not. present. The method is 
versatile enough to allow for rectangular fields, 
irregularly-shaped fields and the placement of beam blocks 
or compensators in the field. The use of measured data as a 
data base is abandoned and replaced with data generated from 
first principles by Monte Carlo techniques. 

The Monte Carlo method is used to map the spatial 
distribution of charged particle energy away from a primary 
photon interaction site. This, distribution, called a 
"primary dose spread array", is convolved spatially with the 
kinetic energy released (KER) at all the primary interaction 


sites to yield: ithe {iprimary. wdose’ distribution. ny) a 
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heterogeneous phantom, the dose spread arrays need to _ be 
modified. 

Any current method separating primary from scatter dose 
can be improved with the methodology proposed here. Since 
convolution is a well-understood mathematical technique, 
there are conceptual and practical advantages for 
implementing a convolution framework for the spread of 
energy due to. scattered photons: as well as-~ charged 
particles set-in-motion. 

The proposed method was based on concepts which will 
not become obsolete as improved computational power becomes 
available. The fundamental reason for this is the sound 
physical basis of the Monte Carlo and convolution 
procedures. Another reason is an inherent flexibility to 
compromise between speed and accuracy. Therefore, as 
computers become faster, greater accuracy can be achieved 
without modifying the algorithm substantially. As well, the 
present availability of "array processors" is compatible 


with the array structure of convolution mathematics. 


* Dean (86) has used this approach to find where the scatter 
dose was being absorbed from Cobalt-60 primary photons 


interacting in one location of a water phantom. 
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7.2 Primary Dose Spread Arrays 


7.2.1 Definition Of A Primary Dose Spread Array 
* 

A primary dose spread array is the 3-dimensional 
spatial distribution of energy deposited by electrons and 
positrons which spread from the site of the primary photon 
interactions. It is generated by tracking the motion of 
ehnarged “particles away \itrom- “primary nteraction sites 
occurring within a cubic volume element (voxel). These 
charged particles are followed using the Monte Carlo method 


in a homogeneous phantom consisting of voxels with the same 


size, atomic number, and density. The amount of energy 
deposited at, and in the neighborhood of, the interaction 
voxel is scored. Figure 99 illustrates the procedure 


schematically. 


7.2.2 The Generation Of Dose Spread Arrays 


Using The MOCA Monte Carlo Code 


The dose spread arrays were generated by a "home-made" 
Monte Carlo code called MOCA developed from the program 


Buildup3.for, the Monte Carlo code used to calculate the 


* The "dose spread array" is analogous to the "point spread 


array" used in image processing. 
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Figure 99. 
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[voxe Dimension, 2 
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Dimension, 2 


Primary photons interact in the inter- 
action voxel and the charged particles 
set in motion are followed through the 
phantom to see where they deposit their 
kinetic energy. 
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build-up dose in Section 3 of Volume 1. MOCA takes into 
account the photoelectric effect, Compton effect and pair 
production. MOCA produces 0.511 MeV photons when the 
positron annhilates. Like Buildup3.for, MOCA employs’ the 
approximations of continuous slowing down (without 
generation of secondary "knock-on" electrons) and Gaussian 
lateral scattering for charged particles. Bremsstrahlung 
was not included in MOCA. Repures [00° to (102) ares flow 
charts of the program MOCA. 

Ae photon history begins, in MOCA, by initializing ‘its 
energy, position and direction. The distance the photon 
travels before it interacts is first determined. The photon 
(or the scattered photon created by a Compton interaction) 
is followed until its energy falls below the photon energy 
Cut—oi-t...or---i.t...inieracts) bythe photoelectric or pair 
production effects, or it leaves the phantom. The type of 
pheton.-interactionuis.chosen by finding “partial fractions” 
based on the ~Compton and = Ppadr 7) production attenuation 
coefficients. Charged particles generated by photon 
interactions are placed in a "queue" for further processing 
after the end of the photon history. 

MOCA employs a novel technique to determine the kinetic 


energy acquired by a charged particle following a Compton or 


a pair production interaction. The Compton interaction 
cross-section per unit kinetic energy of the recoil 
electron, do/dT , has a complex dependence on kinetic 
energy, T. The method provides a way of numerically 


inverting the differential cross-section with respect to 


jest S SOS CLG oe SP “e: Bi: 


+n? bale wee peg m2. 


fsish. tay he  arnmaeset at een te 


eo aate 
rae 4’ 


srt 


rer y 


7 seedy! | talk kbgy abe, ‘ohh es Pil 


os ort Prey a 


ate a Ay: 


.o 


t6h oF ope Priogds Tavad | s axoTaqne 4200" 


tint 9foisea@ Bsgoato & YG, betidpms YRTeNs 7 


» A = J cae: —7 


J te " 


ame Tov? a # passe ar : eat ae 


ao 


= a 


rqdtod See bcabadbric:, ‘as 


2 a 


] 
S 


fot 150 Geiaihion s20M.: = 1 
‘eiiitad 44d senette 
siivots sRenmt fo0S “ to.) ete 


es oa ih 


[329 69 "deeaaigna” ys shdesee doin wy 
51.) wae bewtale (ee? faviret snoe | 
GT esainget | AON; ai” ‘peputont 

a «hoe NeTRZOIW ea 10 9 
yd een! Fi aah eeie Yotatd seta, § 
b rSsnt yy ee bon notiteog. 


Ye ed nial bores 269 aS 
a ae fivae ee 


I woleod ahiky 


iain tees ‘Bab vam 


wild 
Tat* rete By fl Aematg o8 top tseanejal 


iysote pit cto att 1S hire aid setta 


Oe cd? ° aa toatatnt gelmiboa iheq 8 ; 

vatom “bitentd dune sot AG PLAee+ eRe Ga 
sirennaoeh stake * Zan 7 “Te\ob: Mattoels 

Yew 8 pasty ree ben een. oT oe YRvens 

| self Mr = 20 rs be Praia stn ees Paitrevn! _ 


514 


Reginning 
Of History 


Reset The Charged Particle 
Counter and Annhilation Flag 


A) 
8) 


Find Out The CT Number 
Of The Volume 


Interpolate The 
Attenuation Coefficient 


Calculate The 


Interaction Distance 


Does The Photon 
Interact In The Voxel ? 


Find The Shortest Distance 
To The Next Wall 


Is The Photon Fnergy Y 
Below The Cut-Off 


? 


Has The Photon Reached 
The Phantom Boundary ? 


Transport The Photon 
To The Voxel Wall 


yf 
Y Is The Photon Fnergy Between 
0.2 MeV and 3.0 MeV ? 
N 


Interpolate The Compton And Pair 
Production Attenuation Coefficients 
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(D) E) 


Figure 100. 


"Partial Fractions" 


Flow chart for photon transport part of 
MOCA. A) and B) are from Figure 101. 
B) can also come from Pagure: LO2ZsjpeC 
D) and E) enter Figure POW. 
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Calculate The 
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Figure 101. 
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Calculate The 
Scattering Angles 
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Of Particle 


ron 


Has The Positron 
Parameters Reen 
Calculated ? 


Cc 
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Has An Annhilation 
Event Occurred ? If 
So, Follow The Second 
Annhilation Photon 


Determine The Emission 
Angles, And Recall The 


Positron Position And 
Assign Photon Fnergy 


F) B) 


Flow chart for the part of MOCA dealing 
with the type of photon interaction. Ca 
D) and E) are from Figure 100, A) andes) 
enter Figure 100. F) enters Figure 102. 
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F) 


Have All Charged Particles Y End Of 
Set In Motion Been Followed ? ; History 
N 


Recall Kinetic Energy, Scattering 
Angles, And Position 


Is The Kinetic Energy 
Less Than The Cut-Off ? 


Deposit The Remaining 
Energy In The Voxel 
Was The Particle 
A Positron ? 
Y 


Determine The Emission 
Angles And Assign The 


Calculate The Kinetic Fnergy 
Lost And The Path Length 


Calculate The Position Of 
The End Of The Step 
Energy To The 

Annhilation Photon 


Has The Charged Particle 
Left The Phantom ? 


Deposit The Fnergy In The Voxel 


Figure 102. Flow chart for the charged particle trans- 
port part of MOCA. F) is from Figure 101. 


B) enters Figure 100. 
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Kinetic energy. 

Figure 103 illustrates the Compton differential 
eross=section “as ay function of kinetic energy up (to the 
maximum kinetic energy, Tmax. The total area under the 
curve is equal to the total Compton cross-section, oO ° 
The curve in Figure 103 is subdivided into N intervals’ such 
that the area of each interval has an equal area G/N 
Each interval represents an equiprobable occurrence of 
interaction. The corresponding kinetic energy of each 
interval, T, (hv), is found. Since each T, (av) is equally 
likely to be the recoil energy of the electron, there can be 
a one-to-one correspondence with a random number. LAs. rs 
accomplished by normalizing the interval kinetic energy , 
T, (hy), to the maximum kinetic energy. 

Figure 104 illustrates the values of T(hv)/Tmax as a 
function of random number, R, for various values of incident 
photon energy, hv(MeV). It can be seen that at higher 
incident photon energies there is a greater probability of 
the generation of high energy recoil electrons. 

The pair production interaction is treated in the same 
manner as the Compton interaction except that the maximum 
kinetic energy, Ty» is now given by Equation 3.1.10. Figure 
105""?4llustrates the’ Cvalues?rotf T(hy )/Ty as? aldunet ton of 
random number for 3.0 MeV and 10.0 MeV. photons. The most 
probable kinetic energy of a charged particle is around half 
the maximum energy. 

The a aatreed particles generated by | photon are taken 


from the "queue" at the end of the photon history. Their 
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Figure 103. ~The Compton differential cross-section as 
a function of kinetic energy. The area 
under the curve, equal to the total inter- 
action cross-section, is divided up into 
equiprobable areas. 
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Figure 104. 


T(hV)/Tmax aS a function of random number, R, 
for various incident photon energies. The 
Compton recoil electron kinetic energies can 
be chosen with equally weighted random numbers 
chosen between 0 and 1. 
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Figure 105. 


/Tp 


T(hv)/T, as a function of random number, Ry 
for various incident photon energies. The 
kinetic energy of pair production charged 
particles can be chosen with equally weighted 
random numbers chosen between Q and l. 
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kinetic energy, position and direction in which they were 
set in motion is recalled. The charged particle transport 
iseaimost 1dentrcal to, Buildups.for. The exception is that 
annhilation photons are produced when a positron reaches the 
end of its range ("annhilation in flight" is ignored). The 
annhilation photons are treated as a scattered photon. 
Appendix 10 contains a listing of the program MOCA. 
MOCA was used instead of EGS because MOCA was much simpler 
than EGS and, therefore, there was much more control over 
the programming of the geometry to describe the dose spread 
array s.may Despite itis Se oy MOCA produces almost. the 
same results as EGS. Figures 106 and 107 are a comparison 
between MOCA and EGS Monte Carlo calculations of the 
depth-dose from Rs, MeV and 5 MeV photon beams, 
respectively. The field size was 10 cm x 10 cm for the MOCA 
calculation “anda fieldea radius: of 5.6 cm Cthe equivalent 
circular radius) for “ay 10) cmax™lOvem. field)s for the EGS 
calculation. Also included in Figure 106 is the measured 
depth-dose for a Cobalt-60 beam. There is good agreement 
between the Monte Carlo programs at both energies and with 


measured data. 
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Figures 108 and 109 show voxel elements of primary dose 
spread arrays for a pencil beam of 15 MV x-rays interacting 


in a homogeneous phantom with cubic voxel dimensions of 
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Figure 106. 
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A comparison of components of the percent 
dose predicted by MOCA and EGS for Cobalt- 
60 photons (1.25 MeV). The measured total 
dose is also shown for comparison. 
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Figure 107. 
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A comparison of components of the percent 
dose predicted by MOCA and EGS for 5 MeV 
photons. 
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Figure 108. Primary dose spread arrays for 15 MV photons 
interacting in an interaction voxel (shown 
with bold boundaries). The voxels are 
cubic with dimensions of 1.0 cm on each 
Side. 
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Figure 109. 
(see Figure 108). 
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1.0 cm. The incident photon fluence spectrum obtained using 
EGS in Section 6 was also used by MOCA. The density of the 
arTraysSi nds © from 0) Preis to ee ca in steps of Owevorse 
The media are assumed to be _ water-like in chemical 
composition. The voxels shown in these figures are in the 
same plane as both the interaction voxel and the _ primary 
pencil beam. The beams have the same cross-sectional area 
aSarthe svoxelSesand* arenomade “ito sinteract™ onlyavin® the 
interaction voxel (bold borders). At least 10° charged 
particles/cm~ were set. in motion by primary photons to 
generate these arrays. To speed up the convolution process 
for rectangular fields, the dose spread arrays are produced 
in a 3-dimensional Cartesian geometry. The dose spread 
arrays are quadrilaterally symmetric about the primary 
photon direction, which is designated as _ the Ak axis. 
Therefore, an element of a dose spread array can be given 
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The dose spread arrays are also symmetric with respect to 


intjemehange.-of the jis iandinAj iaxisyrtherefiore: 


A,(0,2,Ai,05,0k) = A, (02,45 04,04) 
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assumption of local energy deposition at 15 MV is untenable. 
The value of the dose spread array at the interaction voxel 
(Tocat Lon 0,0,0) i nea 1.0g/cem? phantom with a voxel size of 
Pom is 0.329 which indicates 32.5 % of the energy released 
at the interaction site is deposited there. The value at 


the same location for a 0.2¢/cm> 


phantom with a voxel size 
of foem is fonly Oali2. The maximum ilongitudinal range of 
charged particles set in motion in a 1.0g/cm> phantom is 
apowt 74: to 16 cm ands ing a Oe ea phantom the maximum 
longitudinal range is about 20 cm! In the lateral direction 
charged particles can contribute dose about 1 to 3 cm and 5 
to. 7 em from the interaction site for eoelen: and 0.2¢/cem> 
phantoms, respectively. Figure 110 illustrates the primary 
dose spread arrays for 6 MV x-rays and a Co-60 (1.25 MeV) 
photon beam. The spectrum used to determine the 6 MV 
build-up curves (see Figure 69) was used. The density of 
the phantoms was Ose vem and the voxel dimension was 
1.0 cm. As the energy decreases, more of the primary energy 
is deposited in the interaction voxel and the longitudinal 
and lateral range of the charged particles is reduced. At 
Zor Mev jin: a phantom of density A eae the primary 
energy is deposited within 1.0 cm from the primary 
interaction site. Therefore, the assumption of local 
primary energy deposition is more reasonable for Co-60 
beams. 

The primary dose spread array values indicate the 


primary energy deposited in each voxel normalized to the 


collision fraction of the total energy released in the 
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Co-60 (1.25 MeV) 
P=0.5g/cms 


Figure 110. Primary dose spread arrays for 6 MV and Co-60 
photons. The voxel dimension is 1.0 cm. 
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interaction voxel: 


Se ost, G1,43 AK) 
Se aera yeh) ms KER(0, 2,0,0,0) 


C725) 
Where Ap(o,%,Ai, Aj, AK). is the primary dose spread array at a 
distance AL Ay, AK from the interaction voxel. 
TeCp, &, Ai, Aj ,AKDPe™rssethe energy deposited “at that site. 
KER(9,%,0,0,0) is the total energy of charged particles set 
imemotion at the, interaction, site that is lost to. electron 
collisions (not bremsstrahlung) in the phantom (87). 

Since the phantom is homogeneous, the interaction and 
dose deposition voxels have the same density. Therefore, 
the array value is also equal to the dose deposited due to 
electron collisions in the neighboring voxels’ per unit 
kinetic energy released per unit mass, Ke (collision KERMA) 
(87), at the interaction voxel. If the dose spread arrays 
are spatially! invariant (Section? 7.6 «will discuss)’ the 
implications of spatial invariance), Equation 7.2.3 can be 
expressed as: 


ie ee _ DOSE(p, 2,1+41,3+A49 ,k+Ak) 
ioe oth ES) a K.(0,£,1,),k) 


Charla) 


The primary dose spread arrays .can .be produced for 
different voxel dimensions and for a variety of water-like 
phantoms with different densities. Figure 111 illustrates 
the dose spread arrays with a voxel and beam dimension of 


3 
5.0 em and a density of 0.2g/cm . The values of the dose 
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Figure 111. Primary dose spread array for 15 MV photons. 
The value of o-% is the same as Figure 108 
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spread array elements of Figure 111 are nearly equal to 
those of Figure 108 (left dose spread array). Therefore, 
the product of the phantom density, o(efem?y, and voxel 
dimension, &(cm), is a fundamental measure of the _ voxel 
isi Zeya rind junit's tot g/om’. This is just an extension of 
O'Connors theorem (71). O'Connors theorem can now be said 
to be independent of the state of electronic equilibrium. 
This is expected because, for a material of fixed atomic 
number, both the stopping power and the angular scattering 
power are directly proportional to the density (88). Thas 
has also been shown experimentally by Young and Kornelson 
(64,80). 

The primary dose spread arrays are stored for various 
values of voxel size expressed in units of radiological 
depth, ep-% (from 0:2 =0.2¢/em- to 1.0g/cm? in step sizes of 
0.2g¢/em-). This allows’ flexibility inthe» choice of voxel 
dimensions for dose computations and hence the spacing of 
eadculation; » points. For example, a choice of a voxel size 
of 0.2g¢/em" can be used either for a medium with op =1.0g/em> 
fore estab lash. es calkcutation. pointse ispaced, at -avgdtstance.ot 
0.2 cm or in a medium with One ene form,ca lcula.t Lone poandas 
spaced 1 cm apart. When the desired spatial resolution and 
density yields. a: value of 0-%: which, has not, been stored, the 


dose spread array values need to be interpolated from the 


stored arrays. 
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7.3. Scatter Dose Spread Arrays 


The Monte Carlo program used to determine the primary 
dose spread arrays is also used to follow scattered photons 
and charged particles produced by the scattered photons. 
The first scattered dose is scored separately from higher 
order multiple scatter. Themidose i-due:ifto { first Tisca tten 
photons, which is deposited "relatively" close to the 
prumamy: interactioniisite ~irsiiplaced ' Mm iva metruncated <fitrst 
scatter (TFS) dose spread array. The first scatter dose 
deposited "relatively" far from the primary interaction site 
is included with multiple scatter ina residual first and 
multiple scatter (RFMS) dose spread aunay Positron 
annhilation photons are treated as if they were multiple 
scattered photons. Since, at megavoltage energies, first 
scattered photons are mainly forward directed, the location 
of the areata ty voxel within the dose spread array has 
been optimized so that more voxels of the TFS dose spread 
array are scored "down-stream" of the primary interaction 
site than "up-stream". 

The bulk of first scatter photons have been separated 
from the multiple scatter photons because of their different 
spatial distribution) and magnitude of Mtcontributiom tol) the 
total dose . Figures o106: ‘and “LO7?sihave ti lusitra ted ithe 
relative ‘contributions oof primary, first A'scaitwer,, and 
multiple scatter dose for 1.25 MeV and 5.0 MeV photons using 
EGS and MOCA Monte Carlo programs. The contribution of 


first scatter photons is almost an order of magnitude lower 
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than that of primary photons. The contribution of multiple 
scatter at this piteld size is less*than one half that: of 
first scatter at depths to 20 cm. Separating most of the 
first scatter from multiple scatter also allows a choice of 
optimal sizes for each type of dose spread array. Having a 
high spatial resolution when little dose is being deposited 
is not an effective use of computation time (89). The size 
of the TFS dose spread array and the voxel dimension of both 
the scatter dose spread arrays and the primary dose spread 
array is selected to compromise speed and accuracy. The 
more elements in the TFS dose spread array, the slower’ the 
eonvolution calculation. The TFS dose spread array voxel 
dimension can be made smaller than the RFMS dose spread 
array, reflecting its proximity and greater importance in 
contributing to the total dose. The TFS dose spread array 
is >produced; and: stored for’? )the “same o-{ values as the 
Seay dose spread array whereas the RFMS dose spread array 
is stored for values of 0-{.lirom 1.0g¢/em- to 5 Og/ em in 
steps of 1.0g/em~. The TFS dose spread array may be 
truncated without compromising too much accuracy because 
most first scatter dose is deposited close to the primary 
interaction site (77). In order to take into account the 
scatter dose arriving at a field boundary from primary 
interactions at the opposite boundary, the total lateral 
size of the RFMS dose spread array must be at least twice 
thet, of “thes largest field in which the dose is to be 


calculated. 


Pigures 112 and 113 iflustrate examples of the TFS and 
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p = 1.0 g/cm? 
p-@ = 1.0 g/cm? 


Truncated first scatter (TFS) dose spread 


Figure 112. 
i array for 15 MV photons. 
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Figure 113. 
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p = 1.0 g/cm? 
p-@ =5.0 g/cm? 


Residual first and multiple scatter (RFMS ) 
dose spread array for 15 MV photons. 
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RFMS dose spread arrays, respectively. The TFS dose spread 
array has a o°% value of Os enn and the RFMS dose spread 
array 0-2 value is BeOaitemes They are displayed as isodose 
curves. The energy deposited by scattered photons has’ been 
normalized to the same value as the primary dose spread 
array, namely to the collision fraction of the total charged 
particle energy released by primary photons. Therefore, the 
dimensionless numbers associated with the isodose lines 
represent the scatter dose per unit primary “collision "KERMA, 
The TFS dose spread array illustrated has a total width of 
tS iemVand4a.etotal Theight of 15eem. SiAny first ‘scatter photon 
energy not deposited within its boundary is included with 


the RFMS dose spread array. 
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7.4 Convolution Dose Calculation In A Homogeneous Phantom 


The dose spread arrays can be viewed as the response 
throughout all the voxels in the phantom to a primary photon 
Fimpulse™ occurring in one voxel. The ‘dose spread arrays 
can, theretore, be wsed as the kernel’ ina convolution 
calculation to produce’a 3-dimensional dose distribution. 
The dese contribution, DOSE CONTRIBUTION (i+Al,j+Aj, k+Ak) tO 
a dose deposition voxel tr JtAj,ktik due to primary 
interactions at a voxel i,j,k in a homogeneous water phantom 


is illustrated in Figure 114 and given by (see Equation 


Tol ay 


DOSE CONTRIBUTION (i+Ai,j+Aj,k+Ak) = 
A(o+2=2,Ai,j,k) K_(o+2=2,4,5,k) 


(754.4) 


Equation 7.4.1 can be used to generate the absolute dose. 
However, usually only the relative dose distribution is 
required. If changes in the beam spectrum (or "hardening" ) 
is (negliglibie, “KeCo-=%,1,j,K) may be. weplaced by the 


relative primary photon fluence: 


a o(4,i,k 
® (1,3,k) = Sa 


Guede Z) 


* The density, p, in the equations of Section 7.4 and in the 
computer program, Volve.for, is taken to be the gravimetric 


density. 
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Figure 114. Illustration of the dose contribution at 
itAi, j+Aj, k+Ak due to primary photon 
interactions at 1,},k from the inter- 
action point of view. 
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Table 20 Parameters 


Energy, (hvo), (MeV) 


Fluence Weight, FE. 


(u,,/0),  (cm’/g) 


OF G7 Za 6.61 


Lt 11 ala 


0.0857 0.0439 0.0267 


0.0326 0.0241 0.0176 


For Determining ji For 15 MV X-Rays 
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The relative fluence, OCCisj,K),. is) the’ primary “photon 
fluence at the interaction voxel normalized to the incident 
primary fluence at-the central axis. For a parallel ‘beam, 
@'(i,j,k) need only take into account the primary photon 
attenuation: 


-jik 
® G57, K= © GeI,0)e 


CaS) 


Where wu is the KERMA-weighted average attenuation 


coefficient and is calculated using: 


uk : 
Pres S07) (hy) 


u 
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The summation is over all the spectral components (9), where 


Fis the relative weight of the spectral fluence component. 


The energy of the primary beam is (hy ye Uy is the 


attenuation coefficient of the spectral component. Table 20 
Summarizes the parameters for the spectral components. Due 
to hardening, ti varies slightly as a function of depth. A 
constant 0703 ! Pinas used. If the dose at very large 
depths is required, then the effective attenuation 
coefficient should be made a variable as a function of depth 
(i.e. i CZ) oF. Figure 115 illustrates the value of Has a 
function of depth. The value of i agrees well with the 15 
MV effective attenuation coefficients found from zero-area 


TMR data and transmission data (see Figure 116) (90-94). If 
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Figure 115. The variation of the effective attenuation 
coefficient -as a..fumetion of depth. .The 
parameters used to produce this graph are 
found in Table 20. 
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Figure 116. The variation of the measured effective 
attenuation coefficients as a function of 
nominal beam energy. This graph was pro- 
duced from zero-area TMR and narrow beam 
attenuation experiments cited in the 
literature. 
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the beam is diverging, an inverse square reduction of the 
primary fluence must also be included: 


29 Z 
ial tubal) 


OoCisd,.h)a=1 Cerin 0).e D 


CPA 59 
Where, D, is the distance from the source to the interaction 
Site along the primary ray. O"(in i. Kk) Can abso account tor 
the “external contour of ~the “patient “(in “which. .case the 
factor, SSD, in the argument of the exponential in Equation 
7.4.5 would not be a constant) and for beam modifying 
devices such as shielding that alters the relative primary 
fluence, but does not alter the dose spread arrays. 


Including the dose spread array, Equation 7.4.1 now becomes: 


DOSE CONTRIBUTION (i+Ai, j+Aj,k+Ak) « 9° (i,j, kK)A(p-2=2,Ai, Aj, Ak) 
(7.4.6) 


The convolution of the dose spread arrays with the 
relative fluence can proceed in two different ways, called 
the "interaction point of view" and the "dose deposition 
point of view". In general the calculation may be done by 
either method, but there are circumstances’ where one 
approach is more efficient. 

In the interaction point of view, the macroscopic beam 
is composed of a set of contiguous pencil beams each of 
which is followed through the phantom to see where _ the 
primary interactions Occur’. The dose spread array 
contributions throughout the phantom are summed for all such 


pencil beams. The dose distribution in the interaction 
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point of view is: 


DOSE (i+Ai,j+Aj,k+Ak) = 2 2 [2 O' (i,j,k) A(or2=2,A1,Aj,Ak)] 
ay 


"es 
k 

Che decay) 
Where the i,j,k Summation is over all interaction voxels; 
i.e. those voxels inside the phantom and the field 
boundary. Depending on the region of interest, the dose 
outside the beam boundary can also be computed. 

Figure 114 illustrates the calculation of dose using 
the interaction point of view. The relative fluence is 
calculated for each point along a pencil beam (inside square 
brackets in Equation 7.4.7) by following the primary pencil 
beams through the phantom. When a beam blocking device such 
as a compensator is introduced in the beam, some of the 
primary pencil beams are affected. The interaction point of 
view can determine the effect of the changed primary pencil 
beams on the dose throughout the phantom without 
recalculating the entire beam. Appendix “ll “lists ~a 
convolution program, Volve.for, which calculates the dose 
using the interaction point of view. 

The interaction point of view gives the dose throughout 
aearnegion ssof ginterest.) +-Often,only the sdose at asfew,voxels 
is desired. In this case, the dose in these select voxels 
is calculated using the dose deposition point of view. 

In a homogeneous, unbounded phantom there is a 
geometrical reciprocity between the primary interaction and 
dose deposition vonelea Even though the dose spread arrays 


were produced to describe the transport and absorption of 
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dose throughout a phantom due to primary interactions at a 
voxel, they also describe the transport to and absorption at 
a dose deposition voxel due to an equal magnitude of primary 
interactions throughout the phantom. This is essentially 
the same principle as the source-target reciprocity used in 
health physics to calculate the organ dose due to internal 
isotopes (95). 

Figure 117 illustrates the use of the dose spread 
arrays from the dose deposition point of view. The dose 
Spread arrays now represent the dose deposited in the dose 
deposition voxel normalized to the collision KERMA produced 
in the interaction voxels. The dose deposition point of 
view calculation sums the dose contribution at the dose 
deposition voxel due to primary interactions throughout’ the 


phantom. The dose at a point 1,J,K is given by: 


DOSE(1,J,K) « £2 E ®*(I-Ai,J-Aj,K-Ak)A(0-2=2, Ai, Aj, Ak) 
Ai Ajak 
(7.4.8) 


~if the» dose at only a few voxcls ™rsprequired, the dose 
deposition spoints of (view is far@moremeefiicient. Lor 
example, the dose deposition point of view should be used 
when determining the dose along the central axis or along 
transverse profiles. The dose deposition point of view is 


mathematically equivalent to taking the convolution using 


* This reciprocity only rigorously applies to the RFMS dose 


spread array when the phantom is infinite. 
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Illustration of the contribution of dose 


Figure ki/. 
from primary photon interactions at I-Al, 


J-Aj, K-Ak arriving at 1 Wcusingsthe 
dose deposition point of view. 
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Pie Serial  prediicimc90)).. “A disadvantage .to calculating a 
large dose distribution using the dose distribution point of 
view is that a large array for the relative fluence must be 
pre-calculated and stored before Equation 7.4.8 can be used. 

Equations 7.4.6, 7.4.7 and 7.4.8 do not specify which 
type of dose spread array is being convolved. The same 
geeneral equations apply to, all the «dose {spread arrays, 
except the number of terms in the summation and the value of 
the relative fluence will depend on the voxel dimension of 
the dose spread array. 

The convolution method does not take into account’ the 
dose due to contamination. The measured contamination dose 
was added to the calculated dose. In Volume 1 it was’ found 
empirically that the magnitude of the contamination dose 
depends on five parameters; the 2-dimensional position from 
the central axis, the radiological depth in the phantom, and 
the width; a, and length, b, of the field. The amount of 


contamination is given by: 


~2x" ja? ~2y*/b* 
Cle, ¥,z452,b) §2/G(O-z) cYa-b e e 


(7.4.9) 
The amount of contamination as a function of radiologicei 
depth was found from Figure 63 inside the field and from 
Figure 31 outside the field. The mean density, along the 
Hrimarye «ray, strom the <surftace ~of) thes phantom to tine 
calculation point is 0. These graphs have been included in 
the form of a look-up table, G(p-z). The dependence of 


contamination on field width was found to be linear. The 
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constant of proportionality, oc, was 1.0-2/cem for a 15 MV 
photon beam. The effective field size is assumed to be the 
square root sof sienes product, of the, field dimensions. 
Equation 7.4.9 assumes a Gaussian dependence Of the 
contamination on the distance from the the central axis. 
Figure 118 illustrates the measured and calculated 
tissue-maximum-ratio (TMR) for a 15 MV beam as a function of 
depth in a homogeneous water phantom along the beam central 
axis: tor field sizes of-6 cm xaGrcem and 20 cm x 20 cm. The 
calculated TMR curve was obtained by assuming the beam was 
parallel and had no inverse-square reduction of the fluence 
with depth, as required for comparison with measured TMR 
values. The agreement between the measured and calculated 
dose is better than 1 % beyond a depth of aya and 2 -wittirin 
10 %@ in the build-up region. Figure 119 illustrates the 
percent depth-dose curves for a 6 MV x-ray beam for’ various 
field sizes. The agreement is within 1 % beyond care for 


10°.cm x 10_cm and 20° cm x 20 emefields and within & 2% for 4 


4 cm x 4 cm field. 

Figure 120 illustrates dose profiles at ae for. fiela 
Sizes sof: 10) em x lO cm and SO7cm x 30 7em) for a lo MV beam: 
The fluence profile at the phantom surface was assumed to be 
uniform inside the beam boundary and zero outside the beam 
boundary. This assumption appears to be adequate for _ the 
10 em x 10 cm field, but does not account for the "horns" in 
the larger‘’fieddwprofilésactThes fall-offs ainsadose. near the 


beam boundaries due to lateral electronic disequilibrium, 


not geometrical penumbra, is accounted for. Figure 121 
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Figure 118. Measured and calculated TMR's for a 15 MV 
beam as a function of depth along the 
central axis in a homogeneous water phantom. 
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Figure 119. Measured and calculated percent depth-dose 
data for a 6 MV beam as a function of 
depth along the central axis in a homogeneous 


water phantom. 
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Figure 120. Measured and calculated dose profiles at 
dnax for a 15 MV beam. 
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Dose profiles at dy4, in homogeneous phan- 
toms with various densities. The dose is 
normalized to the central axis dose ina 
unit density phantom. 
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illustrates the dose profiles in homogeneous phantoms with 
various densities. The dose is normalized to the dose in a 
unit density phantom at the central axis. The smaller the 
density, the smaller the dose inside the field. This= is 
especially pronounced near the field boundary. The dose 
outside the field is greater for smaller densities. This 
effect has been observed by Kornelsen and Young for 10 MV 
x-rays (80). It is due to the enhanced range of charged 
particles in low-density materials. The reduced dose inside 
the field is mainly due to the charged particles migrating 
inward from outside the field and the increased dose outside 
the field is due to an increased number of charged particles 
"streaming" there from inside the field. 

The method can be used to determine the dose in a 
wedged field. It was assumed that the wedge could be 
described mathematically as a linear (or "ramp") increase in 
the primary fluence from one side of the beam to the other. 
The primary fluence profile is shown in Figure 122 along 
with the tfluence calculated.for the 45° and Be wedges 
(defined by the angle the isodose lines makes with a line 
parallel to the phantom surface at some specified depth) 
supplied by the accelerator manufacturer. The wedge 
material is steel with an effective attenuation coefficient 
of 0.284 emi (97). The closest fluence profile curve to 
the linear profile curve is for a 60° wedge. Figure 123 
illustrates the calculated isodose curves for a 15 MV beam 
with a fluence profile as shown in Figure 122 and the 


O A : 
measured isodose curve for a 60 wedge. A diverging beam 
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RELATIVE FLUENCE 
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The relative primary fluence profile for 
15 MV beam wedges and the fluence profile 
used to obtain the calculated isodose 
curve in Figure 123. 
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Figure 123. The calculated isodose curve for the fluence 
profile shown in Figure 122 and the measured 
isodose curve for a 60° wedge for a 15 MV beam. 
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was used with a field size of 10 cm x 10 cm defined at the 
surface of the phantom which was 100 ecm from the _- source. 
The calculation predicts the characteristic reversal of the 
slope of the isodose lines in the build-up region. The 
agreement of the position of the isodose lines is good and 
might improve by a more closely matched fluence 
distribution. 

The method can also take into account 
irregularly-shaped fields. Figure 124 illustrates’ the 
effect on the dose profile at a depth of 5 cm when a_ shield 
is placed in the 15 MV beam. The shield consisted of a bar 
of cerrobend (lead-tin-bismuth compound ) extending 
completely across the field in a transverse direction to the 
profiles, Die, width of the bar’ was 1.5-cm,j;and, it. shad= «sa 
thickness of 8.5 cm. The primary fluence transmission 
through the cerrobend was calculated to be 3.1 %. Figure 
124. also ilblustrates a calculationwof the profile assuming 
the primary dose is deposited locally. The measured profile 
agrees. best wathy the calculation using non—local primary 
energy deposition. This illustrates that charged particles 
are streaming into the region under the shield. Also 
illustrated in Figure 124 is a calculation of the dose using 
a Clarkson scatter summation technique (98). This involves 
finding the “scatter? by summing. contributions | from the 
unshielded portion ofthe field*and then adding this ‘to the 
primary.dose under »the,.shield, which. is» equal ©to the 
transmission factor (in this case 0.031) multiplied by the 


zero-area TMR. This method predicts about one-half of the 
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Figure 124. 


DISTANCE FROM CENTRAL AXIS (cm] 


The measured dose profile at a depth of 5 cm 
when a shield is placed in the 15 MV beam. 
Tel relative position of the shield sisiin— 
dicated. Also shown are various calculated 
dose profiles for this situation. 
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dose under the shield compared to measured data or the 
convolution method. This is because this method assumes a 
total absence of primary dose under the shield being 
transported from the unshielded portion of the field. This 
method, which has an underlaying assumption that the primary 
energy is deposited locally, agrees well with the 
convolution method in which local primary energy deposition 
has been assumed here temporarily for comparison (this is 
achieved by replacing the calculated primary dose spread 
arrays with an array in which the interaction voxel has a 
value of 1.0 and all other voxels have values of zero). 
Thelfinding in Figure lt2i that some* of*®jithe "penumbra" 
is due to electronic disequilibrium suggests a _ simple 
calculation method for this shield based on measured dose 
distributions. Figure 125 illustrates that the bar shield 
can be represented by three fields. The middle field 
represents the radiation transmission through the shield. 
The fieldsphave: tos be tilted: “soj that Mmeighboring? field 
boundaries remain parallel. The relative weight of the 
shielded field compared to the unshielded field is equal to 
thes transmission factor. Figure 124 also illustrates the 
calculated dose using this calculation method. It provides 
better agreement with the measured dose. The shield is 
acting like a collimator, Lt is -‘shitelding the ftlueneers sbi 
it *#does “not “prevent<"charged™ particles generated ‘in “the 
unshielded region from being transported to the shielded 


region. 
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Figure s125 {6 Bhe iba shield represented by three fields. 
The central axes of the fields are shown 
with dotted lines. The results of the 
calculation with this representation is 
shown in Figure 124. 
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7.5 Extension to Heterogeneous Media 


Charged particle transport through fa! heterogeneous 
medium is much _ more complex. To be rigorous, the primary 
dose spread arrays would have to be generated for each 
heterogenous situation that could be encountered. The 
number of possible combinations is enormous so that) <an 
acceptable approximation is necessary to take advantage Oc 
data stored in the dose spread arrays generated for 
homogeneous phantoms of different densities. 

Figures 95 and 97 clearly illustrated that in a 
low-density medium, charged particles migrate a greater 
distance from the primary interaction site. Tase musi. be 
taken into account when describing the transport of charged 
particles from a unit density medium such as muscle into a 
low-density region such as lung. Figures 126 and 127 
illustrate a MOCA Monte Carlo generation of the primary, 
first scatter and multiple scatter dose components for 6 MV 
x-rays in homogeneous and heterogeneous medium, 
respectively. The natural logarithm of the percentage of 
dose normalized to the maximum dose is plotted as a function 
OL depth. The primary component decays approximately 
exponentially with depth beyond a max in both the 
unit-density and low-density regions, however, the decay 
constant is smaller in the low-density medium. The multiple 
scatter dose has an extended build-up in the unit-density 
medium. it Oe: decay exponentially in the low-density 


medium. Instead, the multiple scatter component first 
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Components of the percent depth dose in a 
homogeneous water phantom predicted by 
MOCA for a 6 MV beam. The measured total 
dose is also shown for comparison. 
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Figure 127. 
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Components of the percent depth dose in a 
heterogeneous phantom (with water chemical 
composition) predicted by MOCA for a 6 MV 
beam. The measured total dose is shown 
for comparison. 
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decreases with depth after the first interface (until 
10-14 cm beneath the surface) and then increases before the 
second interface. The first scatter dose has features of 
both the primary and multiple scatter components. It has a 
build-up in the unit density material similar to the 
multiple scatter dose but has an approximate exponential 
reduction in dose in the low-density material. 

The extended build-up in the scatter dose seen in 
Figures 126 and 127 is reflected by the large longitudinal 
extent of their dose spread arrays. However, the behavior 
of the first scatter dose with depth in the heterogeneous 
medium resembles the primary dose component. 

The extension of O'Connors theorem for charged 
particles set-in-motion suggests that range scaling may be 
used to approximate charged particle dose spread arrays in 
heterogeneous’ phantoms. The dose spread arrays stored for 
various values of o-2 were, therefore, used. The spatial 
resolution of the dose calculation, %, is chosen and fixed 
ati the start, off.theycalculation. As) shown in, Figure (128, 
the average density, 0 , between the interaction and dose 
deposition sites is found. The: tarray value for the 
location, Ai, 4j,Ak in a homogeneous phantom with the same 
6-2 value is found by interpolating between dose spread 
arrays at? fixed o-2 values. Linear interpolation is 
sufficient because the array values vary relatively slowly 
as a function of p-%. 

There are many algorithms for finding the equivalent 


1 > A 
homogeneous density ‘environment between the interaction and 
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Figure 128. 
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Dose Deposition at 
it+Ai, j+Aj, k+Ak 


Determination of the average density be- 
tween the interaction and dose deposition 
voxels. 
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dose deposition sites. The simplest, conceptually, is 
Pray-tracing’ ~ Samples of density are taken at evenly 
spaced intervals between these voxels. Figure 129 


illustrates the procedure schematically. Suppose that the 
number of intervals is N. The number of samples is N+l 
(including the beginning and end of the path). The length 


of-each jinterval in whe Cartesian directions is oi, 4j,6ék 


where: 
6i = Ai/N 
(7.5.1) 
6j = Aj/N 
G2 <S2) 
6k = Ak/N 
C7te5e. 3) 
The average density is given by: 
N-1 
O-= [£,e (i,j,k) ee o(i+ndi,j+j6j,k+ndk) + fy (1+Ai ,j+Aj ,k+Ak) ] /N 
C7e5isao) 
The weights of the first and last voxel are ee and NE 


respectively. These weighting factors are usually set equal 
to 0.5 because the mean position of charged particles 
set—-in-motion is in the middle of the interaction voxel. 
Interpolating to get a dose spread array valid for an 
average density between the interaction and dose deposition 
sites, 6, implies that it has been normalized to the total 


amount of energy released in an interaction voxel of the 
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Figure, 129. Ray tracing is performed by sampling the 
density between the interaction and dose 
deposition voxels. 
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same density. However, the density of the interaction 
Voxel?..0 Gi, 47k): in a heterogeneous phantom, is generally 
not the same as the average density between the interaction 
and dose deposition voxel. The primary dose contribution to 
a dose deposition voxel, due to primary inveractions): #at? | san 
interaction. vdxel of density oe (i,j,k) when the average 


density between these voxels is §, is given byes 


i oe ah i,j K)AD(6-2,Ai,Aj,Ak 
DOSE CONTRIBUTION (i+Ai, j+Aj,k+k) « 07(i,j,k) a Sue Sisco) 


i. 5 5p) 


Pie hac tor Ped.) , kK )i/e putakes mer ntommaccountmethe: edattencnt 
amount of kinetic energy released in the heterogeneous 
interaction voxel compared to the amount set in motion in 
the interpolated homogeneous voxel of density 650). 
Otherwise, the convolution is carried out in the same manner 
as in a homogeneous’ phantom. Therefore, this procedure 
avoids the need to first compute the dose in water and then 
calculate an "inhomogeneity correction" separately. 

If the phantom is homogeneous, the factor DIGiTs jiaky, Pin 
Equation 7.5.5 becomes 1.0 and the equation simplifies to 
Equatron 974.6. Im/Many\ ‘case Wasi quarts on pinion can be 
Simplified to two terms. Oi ds ky paepends. Only vons tie 
interaction voxel and p depends on the average path density. 


Therefore, Equation 7.5.5 becomes: 


DOSE CONTRIBUTION (1+A1,j+Aj,k+Ak) <« $°°(1,j,k)Ap’ (6-2, Ai, Aj, Ak) 


(7-536) 


sist tive hay ts. ei + phe) SEE. .fax0v pet 


sft jo = tsGh re  nevewall ; 


5! paar 
nei : nova. Shame agenajed at ye P td, : ie ce 
on Wi - * : 


on ci af «2 aiem — ant_. on 
ila teen etl betel ume tixdy | sofre08 4 
r at wtenoy sani nop tee. 


(7 


e Mad 


(hea) a at cocina 


3 iy > ia ay 


4 ~ 
2S 
7 { 
4 
r 


hbase a2 silat Va, Lm 2) q noroa 

oF WE © boaseign darance oltenin.. +0 Nune 
dk 9 Aad yea wr ai Resegroo - texov potion. : 
i+ 3G Tene ts | * ) “Bedelearsiat 
ni hua errr a bru Z0uqo. eid oaierat a6 

ng Leaae “i pha cece 6 oni al ; 

a rei 9268 et f “> peti ea Of busca ed? ebtove 
jathaee os et ts ase eee” ns Gyaivae 


"3¥ moinedg ad 11 


> 


mteoto Lisype at! best —egmonad 8.6.7 noltaupa + 
soc ypsups - sse50 ea Ri Be. ‘noltaup® 


Ving wehasgeor (NE, 204 aries ow? - 3 heltitqmie. 


a 


[k7 no stpeiet.) Bae tenov gotionisia3 


-sagoaed #18.% sei temp S1ototedt 


=~ 


i 


(yt) © GULBe Tish tee Pa nUSENMOO 3200 


368 


Where: 
o°°G,554) = ©°(4,j,k)0(4,j,k) 


Pai) 


Ap(6+£,A1,4j ,Ak) 
Ap“ (6-2,Ai,Aj,Ak) = 


p 
(See) 

The TFS dose spread array may be range-scaled using the 
same algorithm. In fact, range Scaling the TFS dose Spread 
array gives the exact correction. This is because lam Tribals etets 
Scatter photon interacting tat. the. dose deposition voxel 
could not have interacted anywhere other than along the path 
between the primary interaction and dose ‘deposition voxels. 
Since the primary dose Spread array and the TFS dose Spread 
array are stored for the same values Ofe to 8' 'the primary and 
TFS dose spread arrays may be combined where they coincide 
Spatially. The TES dose Spread array need only be dealt 
with separately if it extends beyond the border of the 
primary dose spread ain aye. 

Equations UA and PADIS have two impaled 


assumptions; 


1) charged particles only travel on the Straight-line path 


or "ray" between the interaction and dose deposition sites. 


2) each portion of the path between the interaction and dose 
deposition sites has an equal influence 6 the 
transportation process. 

Both of the above assumptions are, strictly-speaking, 


untrue. Charged particles interact ‘and Scatter almost 
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continuously so they can travel on any path between the 
interaction and dose deposition sites with only the 
constraint that the overall distance travelled is less’ than.- 
the maximum possible pathlength: or, “ranger: Charged 
particle scattering between two points is a highly complex 
process which cannot be described exactly by an average 
density. However, the most probable of any path= “rss tine 
direct one between the interaction and dose deposition 
Sites. Therefore, the assumptions should be viewed instead 
as a first-order approximation to the solution of charged 
particle transport... 

The validity of the approximations were tested by 
computing primary dose spread arrays directly using the 


Monte Carlo method for a variety of water-like heterogeneous 


phantoms to Simulate transport across tissue-lung, 
lung-tissue and air gaps. Figure 130 illustrates, 
schematically, the phantom and the? position Jof the 


interaction voxel (the jvoxel with bold borders) in the 
phantoms that were tested. Dose spread arrays generated in 
the heterogeneous phantom from the Monte Carlo method were 
compared to . the. )calculated dose spread arrays’ for 4a 


heterogeneous phantom uSing: 


OO (Ai Aa N= beet ne (Ot REI) 


Dtist ptlst 
G/s5.9) 
Where Meer ot CPE At Ady BK) is the primary and TFS dose 
Spread arrays combined. Tables. A. stOs.dse4 ims Append iaxé 12 


compare the results of the two calculations. The dose 
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Figure 130. 
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Schematic representation of the slab phantoms © 


tested to verify the approximations used in 
determining the dose in heterogeneous phan- 


-toms. The interaction voxel is shown as 


a square. The results are found in Appendix 
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Spread arrays are in close agreement with respect to the 
range and absolute value. The maximum deviation between the 
two values occurs near the heterogeneity interface. A check 
was also made to see if the method chosen to find the dose 
spread arrays in Equation 7.5.9 conserved energy. Equation 
7.5.9 gives the amount of dose/unit KERMA being deposited at 
the dose deposition site, therefore, the sum of energy/unit 


KERMA is given by: 


ENERGY SUM = 5 5 fA 9 (Ad Aqsa) oCisAi 4A) Rak) 


Aiajak Pttst 
C725 510) 
The value obtained from Equation 7.5.10 was compared to the 
result from the Monte’ Caro" [calculations using the 


following: 


; -«n, . MONTE CARLO SUM - CONVOLUTION SUM 
* DEVIATION OF ENERGY SUM = ST TI ORe TO ROSIN) ok 


(iw 5' al 1) 


The comparison is shown at the bottom of Tables A to I in 
Appendix 12. Energy was usually conserved to better than 
5 %, except for Situation G, in which the percent deviation 
was 8 %. The results can be improved by modifying the 
ray-tracing algorithm. Instead of giving f a weight of 0.5 
for all rays, some rays are given an fy weight greater than 
OO. Tables -A@ “to L” “sin Appendix sl2-  Sillustrate.seine 
heterogeneous calculation with oe equal to 0.9 when Ak is 
positive and Ai and Aj are both O and fy equal se) 0.6 when 
Akers posit myevandayimores}) COuUtmnOU both) is) +l Al eother 


f values are still 0.5. Energy is now conserved to better 
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Phan 2) % for vaile values except those in Situation G for 
which energy was conserved to 95.5 % The results would 
probably also _ improve Lf the value of fy was not constant. 
This has not been tested. 

| Multiple scatter photons may have interacted anywhere 
within the phantom, not just between the primary interaction 
and dose deposition voxels. Therefore, the average "global" 
density of the phantom is usede in Equation, 7.0.0 and if ar ane: 
instead of the mean density along the path between the 
interaction and dose deposition voxels. This approximation 
is. further justified: because multiple scattering is a 
minority contributor to the total dose. This approximation 
avoids calculating the mean density between the interaction 
and dose deposition voxels of the RFMS dose spread array 
resulting in a saving of computation time. 

Figures 131 and 132 illustrate the IMR correction 
factor for a 15 MV beam along the centralcaxis fora 
heterogeneous phantom consisting of horizontal slabs Or 
low-density and unit density material for field sizes of 
5 cme x 0 cm and 10V%emex, 0 ‘cm. The experimental 
measurements were the same as those obtained in Section 6. 
The correction factor is less than On cork at 
5 em x 5 cm which indicates that the dose to cork is less 
than the dose in the homogeneous phantom even though the 
primary photon fluence is greater. This was shown (in 
Section 6) to be due to a loss of lateral electronic 


equilibrium because the distance from the central axis to 
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Frgure 131. the experimental and measured TMR correction 
factor for a 15 MV beam. The experimental 
set-up is snown, in Figure 85. “The freltd 
size srs, Oo. cm x9 CM. 
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Figure 132. The experimental and measured TMR correction 
factor for a 15 MV beam. The experimental 
set-up is shown in Figure 85. The field 
Sik7eniselOnecmmx GLO ¢cms 


oltar RAT > 
\ , pede 2 
¥ ae nat © 
a 
owiesesht 
i i ao " 
4 p , 4. H 


sy & we EMI 


Sree as i: Fe OUT 


: s a 


“avg OF = @ 


ST agro 


th stab er Sig 


LGa qxinged ats 


we OF 


' 
; qU~i pe 


eT et 


a 
a 


ze oat 
orna’ 


asi¢ 


anugit 


the field boundary is smaller than the lateral range of 
charged particles. The dose results of some existing 
methods are shown for comparison in Figures 131 and 132. 
The existing methods assume electronic equilibrium and 
poorly predict the dose in this situation. The convolution 
calculation predicts not only the correct trend but also 
predicts the measured correction factor to within 2 &%. 
Lateral equilibrium is established at a field. size of 
lO -cmexutO/cm,/ and | for this). Situation,’ the convolution 
calculation and the equivalent TAR methods both predict the 
dose adequately. 

The choice of weighting of the interaction voxel, fo? 
ine “the, ray-tracing algorithm’ -altered. the dose’ in this 
heterogeneous situation by only a few percent. Therefore, 


at least for the slab geometries tested, ray-tracing seems 


to be adequate for determining the average density 
"environment" between the interaction and dose deposition 
voxels. 


Many dose calculation algorithms, such as the Batho and 
TAR ratio methods, only correct the dose at the central axis 
for tissue heterogeneities. Figure 133 illustrates the dose 
profile at a depth of 9.5 cm inside a 0.30g/cem> region of a 
heterogeneous phantom (see inset). The field size is 
Z20ecm X..20 em. Atuhthesl centrale a xis 9 the Vadose sin this 
heterogeneous phantom is about 5 % greater than the dose in 
a homogeneous phantom at the same depth. However, the dose 
in the heterogeneous phantom decreases faster than the dose 


in a homogeneous phantom at greater distances from the 
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Figure 133. The calculated homogeneous and heterogeneous 


TMR profile at 9.5 cm depth (A = 4.5 cm) in 
a 20 cm x 20 cm field. 
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central axis. Just inside the field boundary, the dose in 
the heterogeneous phantom is less than the dose in a 
homogeneous phantom. Outside the field boundary, the dose 
in the heterogeneous phantom is once again greater than in 
the homogeneous' phantom. The behavior near the _ field 
boundary, .in the low-density material, can be attributed to 
charged particles streaming from inside the field (reducing 
the dose there) to outside the field (increasing the dose 
there). Therefore, an inhomogeneity correction facroxr 
obtained at the central axis cannot be applied throughout 
the field. Kornelsen and Young (80) have measured similar 
changes in 10 MV dose profiles in phantoms containing 
low-density regions. 

Traditionally, the effects of a wedged field incident 
on a heterogeneous phantom on a dose distribution has been 
performed in two steps. Tbeo doses distribution. es.) first 
calculated for a wedged field incident on a homogeneous 
phantom. "? The correetion. factor sis fthen obtained “for a 
heterogeneous phantom irradiated by an open field. The 
inhomogeneity correction factors are then multiplied by the 
dose distribution for the wedged field to obtain the dose 
distribution. This approach is non-rigorous and has never 
been tested for its validity. The convolution method is 
capable of determining the wedged field dose distribution 
for a heterogeneous phantom in one step. 

Figune 164 > iitmsiratest the tiritheres "is" > good agreement 
between separating the tasks of determining the wedge dose 


distribution and the heterogeneous correction and 
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Figure 134. A comparison of separating the task of deter— 


mining the wedge dose distribution and cal- 
culating an inhomogeneity correction factor 
and performing the calculation in one step. 
Thesresult in erther case is the percent 
dose profile. The field size is 10 cm x 

10 cm. The heterogeneous phantom and wedge 
are shown schematically. 
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calculating the wedge and heterogeneous distribution 


directly in one step. 
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7.6 The Spatial Invariance Of The Dose Spread Arrays 


The dose spread arrays are produced using a very large 
phantom. The primary dose spread array is finite in extent 
so the phantom need only have been larger than the maximum 
range of charged particles set-in-motion. The TFS dose 
spread array is completely independent of the size of the 
phantom used in the simulation. The size of the phantom 
will affect the RFMS dose spread array only. The ‘use: -of 
this dose spread array near the edge of the phantom will 
introduce error. However, since this dose spread array 
contriroutes, (a) small fraction to the, total dose ,. this: error 
due to the extent of the phantom, will be insignificant. 

The dose spread arrays are not completely spatially 
invariant in homogeneous’ phantoms’ for two other reasons. 
Beam hardening is due to the preferential removal of low 
energy photons from the beam. It is possible that the beam 
will have a "hardened" spectrum near the central axis, due 
to a greater pathlength of the primary beam through the 
field’ flattening filter, as. compared to -near theivybeam 
boundary. In addition, as the beam penetrates the phantom, 
it-contains a greater fraction of “higher energy photons. 
Beam hardening effects were not included in the results 
described. If it is necessary to include hardening effects 
for some beams, the convolution can proceed using different 
dose spread arrays for different locations in the phantom. 

Divergent baths have three geometrical characteristics 


which distinguish them from parallel beams. There is an 
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inverse square reduction of the primary photons and a linear 
increase in each field dimension as a function of depth. 
Another characteristic may have to be taken into account for 
this method. The dose spread arrays are not invariant for 
diverging beams. Figure 135 illustrates that a dose spread 
array §should =) be “tilted” at an angle from the central axis 
equal to that of the primary beam. For example, the angle 
of tilt is about 8.5 degrees at a distance Ofte LOOmeme trom 
the source at the edge of a 30 cm x 30 cm field. To” {be 
rigorous, the subscripts of the dose spread array at 


it Ai ,jTAj,ktAk shouldbe transformed to itdAi',j+Aj',k+Ak', 


where: 
ieee (SSD+k) Ak+r Ar 
D 
(i021) 
ee eS AiD ae vA 
SSD+k 
(20a) 
Aj AjD - j4k 
SSD+k 
(Oe) 
/ ; 
eo eee and Araeyiigee Nicos elt hi’ i) one ak lemane 
not integers, interpolation is necessary. This hasenot 


proven to be an important effect when calculating s1MR Sma 
the simple phantoms investigated so far. The transformation 
changes the value of the TMR at the central axis by less 
than 1 % up to a thickness of 20 cm in a 20 cm X 20 cm field 


at a distance of 100 cm from the “source. However, beam 
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Figure 135. 


The dose spread array should be tilted to 
Simulate the dose deposition from a divergent 
primary pencil beam: interacting ‘ait 557k. 
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divergence will be an important effect when calculating 


percentage depth-dose for large fields or may be important 


when calculating dose im phantoms with complex 


heterogeneities. 
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7.7 Comparison With Other Methods and Potential Improvements 


The convolution method presented has features similar 
to the equivalent TAR and A-Volume methods. Before the 
Similarities are discussed, the characteristics of the 


method, which are unique, will be summarized. 


1) The method explicitly takes into account the transport of 
charged. particles. ,which.,allows «he. calculation.of.dose, in 
the build-up region and in lateral disequilibrium situations 
such as near beam boundaries or near low-density 


heterogeneities. 


2) The use of experimental measurements as the data base is 
abandoned. Rather "synthetic measurements" are generated 
using the Monte Carlo method and validated by comparing the 
results of calculations with a limited set of experimental 


measurements. 


3. us Primary »fluence,.is »completely 4 separated... trom dose 


deposition. 


4) The primary dose, local first scatter dose and multiple 
scatter dose can be calculated separately. This allows 
optimization of the calculation resolution for each 
component to find the best compromise between speed and 


accuracy. 


5) O'Connors theorem is explicitly used because the dose 


spread arrays are stored for a range of values of phantom 
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density times voxel dimension. 


6) The same algorithm is used to calculate all the dose 


components. 


Figure 136 illustrates a comparison between various 
dose calculation methods. The Monte Carlo method is the 
only method with the capabilities of the convolution method. 
However, it is yet far too slow for routine treatment 
planning. As computational power becomes less expensive and 
more available, this may change. Nevertheless, 
computational hardware improvements will also increase _ the 
speed of convolution calculations so that the convolution 
method will always remain faster than the Monte Carlo 
method. The convolution method retains the essential 
features of the Monte Carlo method for this application 
without the slow calculation time. 

A procedure common to all other pencil beam and 
pixel-by-pixel calculation methods is" ray-tracing". The 
equivalent TAR method ray-traces to find the density 
weighted with respect to the amount of scatter originating 
at a location. The A-Volume method ray traces to find the 
contribution of first (and some second) scatter dose between 
thescalculation point and all other points, in} the: phantom 
where the primary beam interacts. The dose calculation for 
the first scatter dose contribution in a homogeneous phantom 
can be shown to be equivalent to a convolution calculation. 
This may be shown formally using a continuous 3-dimensional 


integral representation of first scatter fluence. The first 
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Algorithms For Making Dose Corrections In Heterogeneous Phantoms 


Algorithm Pathlength Field Size Position Of Shape Of Flectronic 
Structure Structure Fquilibriunm 

Linear 

Attenuation Yes No No No No 

Coefficient 

Ratio Of 

TAR's Yes Yes No No No 

Effective SSD Yes Yes No No No 

Isodose Shift Yes Yes No No No 

Batho 

(Power Law) Yes Yes Yes No No 

Equivalent TAR Yes Yes Yes Yes No 

Delta-Volume Yes Yes Yes Yes No 

Monte Carlo Yes Yes Yes Yes Yes 

Convolution Yes Yes Yes Yes Yes 


Figure 136. A comparison of the convolution method with 
existing dose calculating methods. 
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> 
Sseatter photon fluence, at- a point r in a homogeneous 
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phantom, due to primary interactions at a DOLD tere Se even 
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Where dp (Tr) is the primary fluence at point r and 0 (2) is 
the density of the primary interaction site. do( ¥-r~)/do is 
the Compton scattering cross-section in units of enaee 

Pal | is the attenuation of the first scatter dose 
between the interaction and dose deposition sites where the 
vector r-r- is directed along the path between the 
interaction and dose deposition sites. Boyer (39) has, shown 
formally that Boauation 7.7.1 is a convolution integral... in 


a formalism used in continuous’ convolution theory. -(96:) 


Equation 7.7.1 becomes: 


bp = | hir)g(t-r-)dV~ 
3D CRE) 


Where: 


h(r*) = op(t)p(r’) 


CS) 


(7.7% <4) 


Equation 7.7.2 is just the 3-dimensional continuous analog 
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Of Equation  ?s4s8. 

In general, the speed of pixel-based ca Veulat ron 
methods is much slower than simpler methods. The time 
depends on the number and type of operations specified by 
the algorithm and the hardware used to run the program. 
Currently, the A-Volume method requires aé_e significantly 
longer time than the convolution method when it is running 
on the same computer (VAX 11/780), presumably under the same 
operating conditions. The A-Volume requires 4.7 hours for a 
PGsxVlG'x°re dose matrix*and 300 *thours for’ a®32< er S32 Sexe 232 
matrix CHM Thée= convoluvien method, by + Céomparison’; 
presently requires 50 min and 8.2 hours, respectively, for 
the same 3-dimensional distributions. The calculation times 
for both of these methods can be improved. The slowest part 
of the pixel-by-pixel methods is ray-tracing. Ray-tracing 
in the A-Volume method is being incorporated into a 
custom-integrated circuit which is expected to increase its 
speed by an order of magnitude (77). The convolution method 
should be particularly well-suited to array processing which 
should"al so "deerease Tethe > "caleula tion by *<an order of 
magnitude. The array processor CAreulat fom -ern tive 
interaction point of view will be examined in some detail. 
The array processor would calculate the following equation 
fora given@Ai ,ige@iktand relative fluence ®! © (i, 77k ye) irom 


Equation *?. S07) 7 
om (1,3,.0A Plo: ’,41,4) ,AK) 
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thes host computerswouldvcalculate © (i,j,k) and also. keep 
track of where’ the primary fluence was interacting in the 
phantom (i.e. the 3-dimensional intersection of the beam 
and the patient). Amat low “chart - of: the host. computer 
calculations is" shown in Figure 137. A flow chart sof - the 
array processor calculation is shown in Figure 138. The 
array processor would first calculate the .mean density 
between the interaction site and all the dose deposition 
Sites. The algorithm then specifies a loop for the 
calculation of Equation Vigaleee? The average density 
calculation could have been brought within the loop and not 
stored as an array. However, some of the average density 
array may be reused for a beam from a different direction 
also} mteracta ngrtait, the same location (i,j,k). Figure 139 
illustrates the idea schematically. The A-Volume method 
should also be able to benefit by storing an average density 
matrix for use with multiple beams. In Surhi, the 
convolution method could employ a custom-integrated circuit 
to calculate the average density array. Both the A-Volume 
and the convolution methods could benefit from an improved 
ray-tracing algorithm. The algorithm presently employed by 
both methods is summarized by Equations 7.5.1 to 7.5.4. A 
more efficient method would be to calculate the density 
along radial vectors. In’ «thes case “of the convohution 
method, in the interaction point of view, the centermor =the 
radial vectors would be the interaction site. In this way, 
the average density, at one point -On). the = ranial vector. 


would use the summation of densities calculated at a point 
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Initialization Of Host 
And Array Processor 


Select The Next Pencil Beam 
E Go To The Next Interaction 

Voxel Along the Pencil Beam 
Calculate The Amount Of 


Primary Fluence Arriving 
At The Interaction Voxel 


Field Completed When 
All The Pencil Beams 
Are Calculated 


Multiply The Fluence By The 
Interaction Voxel Density 
To Get @” (i,j,k) 


Pass On The Interaction 
Voxel Location (i,j,k), Beam 
Direction, And ©”(i,j,k) 

To The Array Processor 


Figure 137. A flow chart of the host computer calculations 
required if an array processor is used to 
perform part of the convolution calculations 
(see Figure 138). 
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Get Interaction Voxel Location 
(i,j,k), Beam Direction, And 
®”(i,j,k) From The Host Computer 


Calculate The Average Densities 
Between The Interaction Voxel 
And All The Dose Deposition 
Voxels Inside The Calculation 
Window And Store The Result In 
An Array O(Ai,4j,dk). If The 
Dose Deposition Voxel Is Outside 
The Window Place A Flag In The 
Array Element 


Go To The Next Dose 
Deposition Voxel 
At End 


Retrieve The Average Density From 
The Stored Array O(Ai,4j,Ak) 
Interpolate The Dose Spread 

Array For The Location 4i,Aj,Ak 


Divide The Interpolated Dose 
Spread Array Value By The 
Average Density 


Wait For The Host 
To Supply Information 
For The Next 

Interaction Voxel 


Multiply By Qa »Jd,K) 


Increment The Dose Matrix At 
Location i+Ai,j+Aj,k+Ak 


Figure 138. A flow chart of the array processor calcula- 
tions (see Figure 137). 


52h 


pond 


yea" oft <4 ried ; tC 
j;camotal vteqet of 
seh 7 


axcal apps sretal 


oe 


,} 


oa Be Be ioeABSoNg ¥% 


5x8 


ChE 


Ay $A aBtza% sect 
Sher, LAF, 


Pealgieitees sr.) 
P “enol? 


adT tnemetsal 
16+ aottacol 


sel ‘erty? 


592 


Common Average 
Density Array Values 


Figure 139. The part of an average density array calcu- 


lated for one beam, that intersects with 
the average density array of another beam, 
may be reused by the second beam. 
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with a smaller radius. If the dose spread arrays were cubes 
with 'n' voxels per side, the number of calculations to 
assemble an average density matrix would be proportional to 
=e with this improved algorithm, compared to the presently 
used algorithm which is approximately proportional to pa 

The proponents of the A-Volume method (77) have 
recently adopted a variable grid spacing capability to speed 
ups ther calenlatiion. tiktaasi phil osophicadil yt isimitartetos: ithe 
convolution method in that it calculates the contribution of 


Scatteriin,ateoarser gridy resolutaon strom contributions: ia 


long distance from the dose point. 
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7.8 Dose In A Non-Water-Like Heterogeneous Phantom 


The present method has been used to calculate the dose 
in water-like heterogeneities such as lung which have a 
Similar atomic composition. The method could be extended to 
other media of different atomic composition such as bone. 

It is instructive to first examine the situation in a 
phantom near the junction of two materials with different 
atomic) numbers. Equations 5.2.1 to 5.2.3 can be extended to 
determine the dose near a heterogeneity due to an atomic 
number change. Suppose the atomic number in a_ region of 
thickness, a, of a phantom is Z and on the other side of 
the interface is Zo + The direction of the photon beam is 
from region 1 to region 2. The dose in region 1 will be 


given by: 
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The first term takes into account the production of charged 
particles) in Pregiong ie that arrive ata depth, z.. and the 
second term takes into account the amount of dose arriving 
at the same depth, z, from the second region. 

The charged particle fluence at the interface will be 
almost entirely generated in region 1 (except for a small 
amount of charged particles directed in the backwards 
direction which has not been included in the above 
equation). There will be a discontinuity of dose at the 
interface. The ratio of dose at the interface in region 2 
compared to region 1 is equal to: 

[S/e]z, 
[S/o)z4 
(728:. 3) 
This discontinuity has been observed by several authors 
K63-3. 54100 )s 

The convolution method may be used to determine the 
dose in medium with different atomic numbers. The primary 
dose spread arrays may have to be recalculated for’ these 
types of tissue because charged particle scattering is 
strongly dependent on atomic number (88). The scatter dose 
Spread arrays would not likely be affected because the 
Compton scattering cross-section per electron is 
approximately independent of atomic number (9). The 
density, 0, would now have to be interpreted as the electron 
density relative to water (67). As in the simple model 
described, by™ Equation. 7.8.1 and 7.8.2,  “twor additional 


factors would have to be included when the primary dose is 
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Z is the average effective atomic number’ between the 
anteraction, and dose. deposition, sites..°. ..he first ratio 
takes into account the relative difference in the amount of 
energy released in the real interaction voxel of effective 
atomic: number, Z(i,j3,k), compared to the amount “released in 
a homogeneous’ phantom of effective atomic number Zn The 
second ratio takes into account the amount of energy 
deposited in the real dose deposition voxel of effective 
atomic number ZCit i,3+ j,k &) compared tos the. amount 
deposited in the homogeneous’) phantom of effective atomic 


number Z. 
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8. Discussion and Conclusions Concerning The 


Convolution Method 


The only justification for our concepts 
and system of concepts is that they serve 
to represent the complex of our experiences; 


beyond this they have no legitimacy. 


Albert Einstein 


(The: Meaning, Of Relativity, L921) 
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8.1 Discussion 


The dose spread arrays are "synthetic" macroscopic data 
determined from the microscopic transport of individual 
particles using the Monte Carlo’ method. Thisy information 
cannot be obtained from measurement because a real primary 
photon beam cannot,.be . "forced" »to interact at only one 
region of the material. This requirement is necessary in 
order to describe the transport of secondary particles’ set 
ingpmotion , by. .primarye.photone interactions... Alimthe dose 
Spread arrays were generated with the inclusion of charged 
pavticle,..transport. The primary dose spread array is 
proportional to the energy deposited by charged particles 
set in motion by primary photons. The scatter dose spread 
arrays are proportional to the energy deposited by charged 
particles set insmotion,by  scattereduphotons. sAntanaly tie 
solution to photon-charged particle transport would be 
@ifEicult (inf not rimpossible)retorobtarns (LO) ey ehorhehis 
reason, the Monte Carlo method was used to generate all the 
dose spread arrays. 

The convolution procedure separates the _ photon - 
charged particle transport up into physically meaningful 
components. Primary -photon tinteractions areusepamated@eirom 
charged particle and secondary photon transport and energy 
deposition .~GLOl).« This «divisions, of srprimary, Bhuence | from 
dose deposition is a more appropriate model of radiation 
transport than the modification of measured dose 


distributions in which these two effects are combined. WAS 
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does not undermine the importance of obtaining good 
measurements, since this provides the selection criteria for 


the photon spectrum used to compute the dose spread arrays. 


The spectrum was determined by first finding a 
tentative "pre-spectrum" which is a measured spectrum 
obtained from a research accelerator (48). The agreement 


between the calculated dose and measured dose is a necessary 
condition for the acceptability of the spectrum, but it is 
not a sufficient condition. The spectrum may not be unique. 
Therefore, a better procedure would be to determine the 
actual spectrum directly by measurement. Transmission 
measurements, obtained with "good geometry", have been used 
to obtain spectra (102,103). A problem with spectra 
determined in this way is that small errors in measurements 
leads to Yarge “errors inf’ (the spectral’ determination. 
Therefore, the spectrum should be checked by getting the 
agreement between the measured and calculated dose using 
this spectrum. If itidid not pagree then the spectrum shnourd 
be adjusted to obtain better agreement (although the 
modification must be kept within the experimental 


uncertainty). 


Transmission measurements are relatively easy to 
obtain. Using a pre-spectrum based on these measurements, 
the convolution method could be adopted as the dose 


calculation method for any accelerator type or energy. What 
is first required is the dose spread arrays calculated as a 
function of p°% and Z for various monoenergetic energies 


from 0 to 50 MeV (no accelerator with a nominal energy above 
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590 MV is available). A spectrum would be obtained for 
energy components for which dose spread arrays have been 
calculated. The dose spread arrays for the spectrum are 
obtained using the pre-stored dose spread arrays. The 
primary attenuation coefficient would then be determined 
using Equation 7.4.4. The dose distribution would then be 
calculated for a wide range of field Sizes. The calculated 
dose beyond qa would be compared to the measured dose. If 
the agreement was not satisfactory, then the spectrum would 
be modified. The above procedure should be repeated for the 
new spectrum. Figure 140 is the flow chart of the proposed 
implementation procedure. 

Basing photon dose calculations on the Prrm 
mathematical foundation of convolution may allow additional 
future improvements. In homogeneous phantoms, the dose 


spread arrays are spatially invariant so that the relative 


fluence array and dose spread arrays may be Fourier 
transformed into the Spatial frequency domain. The 
convolution may then _ proceed much more quickly by 


multiplying the transformed arrays (99). The dose can then 
be obtained by taking the inverse Fourier transform of the 
result. For homogeneous water-like media, this procedure is 


described mathematically as: 
DOSE (x4y;,,2) = ct e{o-(i, j,k) JRFLA(O- 222, Ai, Aj , AK}] 
GSea) 


Where F and ee describe symbolically the 3-dimensional 


Fourier and inverse Fourier transform, respectively (96). 
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Measure The Transmission 
Spectrum 


Assemble A Pre-Spectrum 
With Fnergies Restricted 
To Those For Which 

Dose Spread Arrays Have 
Been Calculated 


Calculate The Dose Spread 
Arrays For The Spectrum 


Determine The Primary 
Attenuation Coefficient 


Calculate The Dose 
Distributions And 
Compare To Measurements 


Satisfactory ? 


Yes 


Modify Spectrum 


Store The Dose Spread Arrays 
And The Primary Attenuation 
Coefficient 


Figure 140. The flow chart of the determination of the 
correct dose spread array and primary 
attenuation coefficient based on the trans- 
mission spectrum of a linear accelerator. 
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This procedure is used routinely in image processing 
applications, and is especially appropriate if an array 
processor is available. 

Deconvolution techniques may be useful in optimizing 
the delivery of radiation therapy. If one knows’ the 
attributes of an optimal dose distribution (eg. uniform 
dose to the target volume and minimal dose to surrounding 
tissue), and the beam directions and intensities, the dose 
spread array can be deconvolved from the ideal dose 
distribution to obtain the best primary fluence 
Gist riput ron. This information could be used to design 
modifying devices, such as beam compensators, to approach 
the ideal relative fluence distribution in homogeneous or 


heterogeneous media. 
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8.2. Conclusions 


The EGS Monte Carlo modelling in Section 6 confirmed 
the interpretation in Section 5 that lateral disequilibriun 
is the reason for reduced dose at. the COneral) ax ise wines 
low-density medium for small field sizes. This situation is 
One aspect. “of a more general phenomenon. Lateral 
disequilbrium never exists near the field boundary penumbra. 
It was shown that the penumbra region for high energy linear 
accelerator beams is due in part to the lateral Cransport .of 


charged particles. When the field size is small and the 


phantom has a low-density, the penumbra extends to the 
Central “axis. In order to describe the disequilibrium 
effects, the transport “of charged particles have to be 


included into dose calculations. 

A dose calculation method based on convolution was 
introduced in Section 7 which explicitly took into account 
charged particle transport. The method was general enough 
to include the transport of scattered photons. The method 
relied on the production of "dose Spread arrays" to describe 
the transport of secondary particles from a primary photon 
interaction voxel to neighboring dose deposition voxels. 
The values of the elements of the dose Spread arrays 
represent the dose deposited in dose deposition voxels per 
UATt Scolltision. KERMA “at a primary interaction voxel. The 
dose spread arrays’ were obtained from Monte Carlo 
Simulations. The primary dose Spread array accounts for the 


primary dose, the truncated first scatter (TFS) dose spread 
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array accounts for the first scatter dose deposition close 
to the’ primary interaction® site,’ and’ the’ residual first’ and 
multiple scatter (RFMS) dose spread array accounts for the 
rest of the first scatter dose and all the multiple scatter 
dose. The primary dose spread arrays generated in 
homogeneous water-like phantoms of density, 9 , with voxel 
dimensions, & , are equivalent to those generated using a 
phantom density,o', with dimensions, ges. provided; 
OL = p%- 2. The scaling also applies to the scatter dose 
Spread arrays. This indicates that O'Connors theorem 
applies to the primary dose even when there is electronic 
disequilibrium, as well as scatter dose, in a homogeneous 
phantom. 

The dose distribution in 3-dimensions is obtained in a 
homogeneous phantom by convolving the dose spread array 
kernel, applicable to the density of the phantom and desired 
voxel resolution, with the relative fluence. There is a 
reciprocity between the interaction and dose deposition 
voxels. The dose spread arrays describe both the transport 
of dose from an interaction voxel and the arrival of dose at 
a dose deposition voxel. This allows the convolution to 
proceed Sineetwoteways Scalledyethe yr rnareraction Said dose 
deposition points of view. The method accurately predicts 
the dose in homogeneous phantoms in situations of electronic 
disequilibrium including the build-up region and near beam 
boundaries. The method was shown to provide aeeunified 
treatment (of dose calculations in that it could account stor 


the changes in dose due to fluence modifying devices such as 
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beam shields and wedges. 
The dose spread arrays, generated for homogeneous 


phantoms at various values of p-£, are used to describe the 


transport of dose in heterogeneous’ phantoms. The average 
density, 0 », between the interaction and dose deposition 
voxels is determined. The primary and TFS dose_ spread 


arrays are interpolated to obtain one appropriate for this 
average density. The convolution of these arrays in a 
heterogeneous phantom contains a factor to account for the 
difference in the amount of interactions between the 
heterogeneous interaction voxel and the homogeneous 
interaction voxel with a density equal to the average path 
density. The RFMS dose spread array in heterogeneous 
phantoms is scaled using the average phantom density instead 
of the average path density. The method provided a good 
prediction of the dose in and near heterogeneous regions. 
Although the method is presently only. applicable to 
dose absorption in materials which are water-like in atomic 
composition, it could be extended to include non-water-like 
materials. The convolution framework of the method lends 
itself to Fourier and deconvolution techniques that could 


improve the calculation speed and optimization of dose. 
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Appendix 10. Listing of the MOCA Monte Carlo Code 


Program MOCA 


c MOnte-CArlo program to generate dose spread arrays. 


¢ The variables are domumented as they occur in the program. 


integer 1,i1,nerg, enn(36), intest 

Teal r,zs,cal,sal,cbeta, sbeta, total 

real stop(179), tha(179), ep, cay 

real eo, e, doset(25, 25, 25), pl, gaus(100),g 
Teal z,max, ang, en(36), dose(5, 25, 25, 25) 


integer is Jsm,p, kay, quad(4), chargpart, inquad(4) 


Teal xnrg(100),mu(100, 10), sig(100,; 10) 
Teal pi(100, 10), que(10), shiftx, shifty 


integer k.nrgint, phantom(25, 25, 25), multi 
integer charge, ir, vint, annhilflag 
integer order, kinorder(25) 

real nrgx,xs,ys,diml, dimlz 

Teal lx, ly, 1z, short, attentot, dist, attensig 
real attenpi, fractionl, fraction2 

real nrg(35), fac, nrgmax, fraci, frac2 

Teal tnrgil, tnrg2. compton(35, 50), tnrg, arg 
Teal kinrg(25), kincal(25), kinsal(25) 
real kinsbeta(25), kincbeta(25), kinx(25) 
Teal kiny(25), kinz(25), kintype(25), rr. 


real pairprod(35,50),v,ro,ri,r2, test(4), arr (25) 


real photox, photoy, photoz 

real oldcal, oldsal, oldcbeta, oldsbeta 
real s,th,arg2 

integer n,numspect 


Teal sum, sqsum,mean(25), stdev(25), sterr(25), maxmean 


real maxsterr, sumpart(4, 25) 
Teal nearx, neary, nearz, near, maxpl, nearo 


integer nearinti, neari, nearinty, neary, nearintk, neark 


integer history, mark, boundary(7) 
integer Oldi, oldy, oldk, smalli, smally, smallk 


common XS) YS) 2S, i, yok, cheta, sbeta, cal, sal, diml,dimlz, short 


c Gaus. dat stores a table of x’s as a function of erf(x). S 


c page 163 and 164 in Volume 1 and Appendix 7. 


c Doseé. dat is a listing of the dose spread arrays. 


c Compton. dat is a look-up table giving the elecron kinetic energies 


c of equiprobable interactions. See Section 7.2. 


c Pairprod.dat is a look-up table giving the electron or positron 
c kinetic energies of equiprobable interactions. See Section 7. 2. 


c PholSnew. dat contains the 15 MV photon spectrum. See Sections 6. 


c and 7.2. 


c Atten. dat contains the attenuation data, scattering and stopping 


Cc powers as a function of energy (from references 9 and 8&8). 


c Debug. dat can contain information for debugging. 


open(unit=1, status=’old’, file=’gaus. dat’) 
open(unit=4, status=’new’, file=’doseé. dat’) 
open(unit=6, status=’old’, file=’compton. dat’) 


412 


& % nck. a 208 eh Gal 


| 


<ONe TIANTEM MOTTAIS 2) .gt tus Ia et ve8 netene? 2 Pics ee 
: i *. “< pie ra Ms reteran sidedorgivps 46a. - aay, 


7 =. Nie 
Otel : I Te) 
UJ i 
“| 1 


Dre) stot Poot sme es, hea 
9 


éun tte tne a2 ‘at ae hons 
#*OoT4s mvt ae ‘ena 4 mob ba asid eres 
a ae 
3 ps chen uttak oe 
ones let «fos ea at ens 


ies? et tiqesa” ‘Leer 
ae 7 


; a bine ei Tle? 
ot see Gam vl « 7 


431 a9; a HM ene saran ts 
jsit Tiivinamneaty 7h 6 grads a 

Pee neon i.vapte segete 
sabaidincbes ih Sab iy atta | re 

ci oie ; Ce el 


$#a077 Seed A 
re te . 


ete rea wedone be 
Slovee seatibae ieebie eabie teoy ” 


‘s2meamte bingy ‘eseng 8S) bee eS), 
Lalenty: :. a 


aan raga COR ee el eee 
sta He? Ripeparast 
“#ED sFR Lo 4 
pimth vies Hight is tee 
ed anny’, af paris r. & eetate teb 2uad > 
rX | hata “AE PRL don COE age > 


were’ Aaitbed Con pn + oi deb evaat 3 


Fiaome 70 = pa iy? cue 53° Hue toot 6 #) ¢54 boverias > 


ar okeee NS eel pyrene sifeadts a ; . 


‘NeLPTA® eee. matsesege posond we 2h. eat Ling sell us si 2 oT 


yee oe prereey ase stab sicepreide edt aeiatrae 16% .nette 
(25 bias 2aans Teyee a ad Fe SIs Eke O80 ‘evguog 2 


sng pedad ae narsaaphat Atetans per ta% poe 2 


. tab agente eee ee en : 
1-jas abe soll Seay ‘euutads dorcnetnes " 
\Yel ao? tae?” eee ‘waudete eet agg 


413 


open(unit=8, status=’old’, file=’pairprod. dat’) 
open(unit=11, status=‘old’, file=’phoiSnew. dat’) 
open(unit=12, status=‘old’, file=’atten. dat’) 
open(unit=13, status=‘’new’, file=’debug. dat’? 


c The following statements read the look-up tables and initializes 


multi=1000 '$ of incident photons 

intest=3 ‘interaction voxel depth indicator 
shiftx=50. 0 !x-position to interaction voxel 
shifty=50. 0 'y-position to interaction voxel 
diml=5. 0 'voxel dimension in x and y direction 
dimlz=5. 0 'voxel dimension in the z-direction 


c Assign CT-number to the phantom voxels. There is a linear relation 
c between density and CT-number (-1000=vacuum, O=unit density) 
c A heterogeneous phantom can be entered 


do 6 k=1,25 
do 6 y=1,25 
do 6 i=1,25 


6 phantom(i, jy, k)=-800 

Cc do 7 k=4, 25 

Cc do 7 j=1,25 

c do 7 i=1,25 

c 7 phantom(1, jy, k)=0 

(3 do 8 k=12,25 

Cc do 8 jy=1,25 

c do 8 i=1,25 

c 8 phantom(i, jy, k)=0 
read(11,345) numspect ‘number of spectral bins 
do 10 i=1,numspect 

read(11,340) xnrg(i),en(i) ‘spectrum weighting 

10 enn(i)=en(i)/100. O#floaty(multi) '# histories/bin 
read(11,350) (que(i),i=1, 10) 'factor used in pair production 
do 12 1=1,100 

12 Tead(1,422) gaus(l1) ‘inverse error function 


‘Nrg’ is the energy, ‘mu’ is the total attenuation coefficient, ‘sig’ is 


the Compton attenuation coefficient, ‘pi’ is the pair production 
attenuation coefficient, ‘tha’ is the scattering power and ‘stop’ is 
the stopping power. 


Nnnnan 


do 14 i=1,24 
14 read(12,570) nrg(i),mu(i,1),sig(i,1),pi(i,1), tha(i), stop(i) 


do 17 i=17,24 
do 16 k=1,5 
16 read (8,550) (pairprod(i, (k-1)#10+m), m=1, 10) 
7/ continue 
do 19 i=1,24 
do 18 k=1,5 


18 read (6,550) (compton(i, (k-1)#10+m), m=1, 10) 
7, continue 
11=95835 'random number seed 
cay=0. 90 'fraction of energy remaining after step 


annhilflag=0 'flags the production of annhilation photons 
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charge=0 ‘number of charged particles to process 
do 200 1=1,numspect ‘do for all spectral bins 
do 180 m=1i, enn(1) ‘do for all photons in the bin 


c Initialize the photon’s energy, coordinate position and angle 
c before transporting. 


xs=ran(il)*#diml+shiftx 'x-position 
ys=ran(il)#diml+shifty 'y-position 

z25=10.0 'z-position 

nrgx=xnrg(l) ‘photon energy of the current bin 
history=historytl ‘photon history counter 

cal=1.0 'cos(alpha); alpha=zenith angle 
sal=0.0 'sin(alpha) 

r=6. 2832#ran(il) 

cbheta=cos(r) 'cos(beta)i beta=azimuth angle 
sbeta=sin(r) 'sin(beta) 

order=1 ‘order of scattering: l=primary-: 


‘and 2=first scattered photon 
c Scattered photons enter at line 20. 
c Test to see if the photon’s energy is below the photon 
c cut-off energy. 


20 if (nrgx .1t. 0.03) go to 50 ‘photon cut-off is 30 keV 
if (order .gt. 4) order=4 'order=4 is multiple scattering 


c Photons which didn’t interact in the last element 
c enter at line 21. Annhilation photons enter here too 


2i continue 


n 


Terminate the primary photon history if it doesn’t interact 
c in the interaction voxel. 


if (ondem aweq. le. andas2s lit. 1070) gol tioslso 
if (order .eq. 1 .and. zs .gt. 15.0 ) go to 180 


c Find out the shortest distance to the next wall along 
c the photon’s path. 


call where 


c Find out the C.T. number of the element and the type of 
c tissue. 


p=1 ‘signifies water 


c Interpolate to find out the appropriate attenuation 
c coefficient. 


nrgint=1 
do 25 while ( nrgx .gt. nrg(nrgint)) 
25 nrgint=nrgintt+i 

if (nrgint .eq. 1) then 
attentot=mu(1,p) 

else 
fac=(nrgx—-nrg (nrgint—-1))/(nrg(nrgint)—nrg(nrgint—1)) 
attentot=mu(nrgint—-1,p)+(mu(nrgint, p)-mu(nrgint—1, p) )#fac 
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end if 
c Find the interaction distance. 


r=ran(ii) 

if (rales (OF0) then 
dist=0.0 

else 
ro=floaty(phantom(i, y, k)+1000)/1000. 0 
dist=-alog(r)/attentot/ro 

end if 


Determine if and if so, where, the photon interacts 

in the element. Rejection sampling is used. Short is the 
distance to the voxel wall in the direction of the photon. 
It was determined in Subroutine Where. 


annan 


if (short .gt. dist) then 


xs=xstdist#sal#cbeta ‘increment the 
yus=ystdist#salesbeta ‘position 
zs=zst+tdist#cal ‘variables 
i=yifix(xs/diml)+1 ‘determine the 
J=sifix(ys/diml)+1 'voxel 
k=yifix(zs/dimlz)+1 ‘location 

c Discard the primary photon and choose another primary photon 

c if it doesn’t interact in the interaction voxel. 


if (order .eq. 1 .and. k .ne. intest) go to 180 


If the photon leaves the phantom, discard it and go on to follow its 
charged particles set in motion. 


an 


if (xs .le. 0.0001 .or. xs .ge. 20. 0#diml) go to 50 
if (ys .le. 0.0001 .or. ys .ge. 20.0#diml) go to 50 
if <Czs) = Fe. (On 000) @and 9 call eelt.2 "6. 0) “gio! to 150 

if (zs .ge. 20. 0#dimlz) go to 50 


Interpolate to get the appropriate Compton and pair 
production attenuation coefficients. The type of 
interaction is determined. ‘’Fractioni’ is the 

fraction of Compton interactions in the total and 
'fraction2’ is the fraction of Compton plus pair 
productions in the total. If the energy is between 

0.2 and 3.0 MeV a Compton interaction must have occured. 


nannannan 


if (nrgx .gt. 0.2 .and. nrgx .1t. 3.0) go to 30 


attensig=sig(nrgint-1, p)+(sig(nrgint, p)-sig(nrgint—1, p)) 
1 #fac 
attenpi=pi(nrgint-1, p)+(pi(nrgint, p)—-pi(nrgint—1i, p))#fac 
fractioni=attensig/attentot 
fractionz=fractionlt+attenpi/attentot 


c Determine which interaction occured. 


rr=ran(il) 
if (rr .1t. fractionl) then 
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c A Compton interaction has occured. The energy 
c of the electron is determined 


30 continue 
nrgmax=(2. O#nrgx/. Sil d¥nrgx/ (2. O#nrgx/. 5ii+i.O) !Compton edge 


nrgint=1 
do 31 while (nrgx .gt. nrg(nrgint)) 
ai - nrgint=nrgint+i 


fracl=(nrgx-nrg(nrgint—-1))/ 
1 (nrg (nrgint)—nrg(nrgint—1)) 
r=ran(il) 
ir=jifix(r#50. O)+1 
frac2=50. O#r—-floaty(ir—-1) 


if (ir .eq. 1) then 
tnrg=(compton(nrgint—1l, 1) 
1 +(compton(nrgint, 1)-compton(nrgint-1,1))#fraci) 
2 #frac2#nrgmax 


else 


tnrgl=compton(nrgint—i, ir-1) . 


1 +(compton(nrgint—-i, ir)-compton(nrgint—1i, ir-1i)) 
2 *frac2 
tnrg2=compton(nrgint, ir-1) 
1 +(compton(nrgint, ir)—-compton(nrgint, ir-1) ) 
2 #frac2 


tnrg=(tnrgi+(tnrg2—-tnrgi)#fraci)*nrgmax 
end if 
c Test to see if the electron is greater than 
c the charged particle cut off energy. If so, 
c then calculate and store the scattering 
c angles and store the kinetic energy and position. 
r=6. 2832#ran(il) ‘azimuth angle in particle coords 
if (tnrg .gt. O. i#ro) then 
charge=charge+l 


kinrg(charge)=tnrg 


arg=(nrgmax-tnrg)#(2. O#nrgx/. Siil+1. O)/tnrg 
g=atan(sqrt(arg)/(1. O+nrgx/. 511)) 'zenith angle 


c Remember the incident photon’s angles in the phantom system. 
oldcal=cal 
oldsal=sal 
oldcbeta=cbeta 
oldsbeta=sbeta 


c Transform from the particle to the phantom coordinate system. 


call geom(r,g.cal, sal, cheta, sbeta) 
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The energy, trigonometric information on the 
angle in the phantom coordinate system, the 
position of the electron and the type of particle 
is stored. 


nanan 


kinorder(charge)=order 
kincal(charge)=cal 
kinsal (charge)=sal 
kincbeta(charge)=cbeta 
kinsbeta(charge)=sbeta 
kinx(charge)=xs 
kiny(charge)=ys 
kinz(charge)=zs 


Recall initial photon’s angles. 


n 


cal=oldcal 
sal=oldsal 
cbeta=oldcbeta 
sbeta=oldsbeta 


n 


If the electron is less than the charged particle 
cut off energy, deposit the energy as dose. 


n 


else : 
To=floaty(phantom(i, jy, k)}+1000)/1000. 0 
dose(order, i» jy, k)=dose (order, i, jy, k)+tnrg/ro 
end if 


c The scattered photon is now dealt with 


order=order+t+l 'the order of scattering is incresed 
c The energy and scattering angles of the scattered 
c photon are determined. 


r=T-3. 1416 ‘scattered photon has opposite azimuth angle 
arg2=1-tnrg#0. S1il/nrgx/(nrgx-tnrg) 
if (arg2 .le. -1.0) then 
g=-3. 1416 
else 
g=acos(arg2) 
end if 
nrgx=nrgx-tnrg ‘determine scattered photon energy 
c Transform from the particle to the phantom coordinate system. 
call geom(r,g,cal, sal, cbeta, sbeta) 
go to 20 ‘scattered photon is followed 


else if (rr .1t. fraction2) then 


c The pair production interaction has occured. First the 
c energy and scattering angles of the electron are determined. 


charge=chargetl 
kintype(charge)=1 ‘signifies an electron 
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Me 


Me 


nrgmax=nrgx-1. 022 ‘available kinetic energy 


nrgint=1 
do 41 while (nrgx .gt. nrg(nrgint)) 
nrgint=nrgint+l 
fracil=(nrgx-nrg(nrgint—-1))/ 
(nrg(nrgint)—-nrg(nrgint—1) ) 
T=ran(il) 
ir=jifix(r#50. O)+1 
frac2=50. O#r-floaty(ir-1) 
if, Gir = eq. 1). then 
tnrg=(pairprod(nrgint-—l, 1) 
+(pairprod(nrgint, 1)-pairprod(nrgint—1, 1)) 
efracl)#frac2#nrgmax 
else 
tnrgi=pairprod(nrgint-1, ir-1) 
+(pairprod(nrgint—1i, ir) 
—-pairprod(nrgint-1, ir-1) )#*#frac2 
tnrg2=pairprod(nrgint, ir-1) 
"  +(pairprod(nrgint, ir) 
—pairprod(nrgint, ir-1) )#frac2 
tnrg=(tnrgit+(tnrg2-tnrgi)#fraci)#nrgmax 
end if ; 
kinrg(charge)=tnrg 


T=6. 2B32#ran(il) 


c This part of the program is done for both the electron 


c and positron. 


is the ratio of total charged particle 


c energy to the energy of the photon. 


43 


nan 


v=(tnrgt+O. 511)/nrgx 
vVint=jifix(v#10. 0)+1 
g=que(vint)#alog(nrgx/. 511)#. 511/nrgx 


Remember the photon’s angles in the phantom system 


oldcal=cal 
oldsal=sal 
oldcbeta=cbeta 
oldsbeta=sbeta 


Transform from the particle to the phantom coordinate system. 


call geom(r.,g,cal,sal,cbeta, sbeta) 


Store the order of scattering, initial charged particle phantom 
angles and position. 


kinorder(charge)=order 
kincal(charge)=cal 
kinsal(charge)=sal 
kincbeta(charge)=cbeta 
kinsbeta(charge)=sbeta 
kinx(charge)=xs 
kiny(charge)=ys 
kinz(charge)=zs 


c Recall the photon’s angles. 


cal=oldcal 
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sal=oldsal 
cheta=oldcbeta 
sbeta=oldsbeta 


c The energy and scattering angles of the positron 
c is determined after those of the electron. 


if (kintype(charge) .eq. 1) then 


charge=charge+li 

kintype(charge)=-1 ‘signifies a positron 
tnrg=nrgmax-tnrg 

kinrg(charge)=tnrg 


r=7Tr-3. 1416 ‘azimuth angle opposite of electron 
go to 43 ‘sent to repeat calculation for positron 
end if 
else 


c A photoelectric interaction has occurred. All of 
c the photon energy is deposited as dose 


ro=floaty(phantom(i, y, k)+1000)/1000. 0 
dose(order, i, jy, k)=dose(order, i, jy, k)+nrgx/tro 


end if 
else 
c An interaction has not occured in the element. 


The photon is transported to the far element 
c wall. 


n 


xs=xs+short#sal*¥cbheta 
ys=ustshort#sal#sbeta 
zs=zst+short*cal 
c A test to see if the photon has left the phantom. 
if (xs .le. 0.0001 .or. xs .ge. 20. O0#diml) go to 50 
if (ys .le. 0.0001 .or. ys .ge. 20. 0#diml) go to 50 
if (zs .le. 0.0001 .and. cal .1t. 0.0) go to 50 
if (zs .ge. 20. O0#dimlz) go to 50 
go to 21 
end if 
c The photon part of the history is completed at line 50. 
50 continue 
c The second annhilation photon is set in motion. 


if (annhilflag .eq. 1) then 


annhilflag=0 
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cal=cos(6. 2832#ri-3. 1416) 'this photon is 
sal=sqrt(i. O-cal¥cal) ‘in the opposite 
cheta=cos(6. 2832#r2-3. 1416) ‘direction of the first 
sbeta=sin(6. 2832#r2-3. 141646) 'annhilation photon 


xs=photox 
ys=photoy 
zs=photoz 


nrgx=0. 511 ‘photon energy assumed to be 0. 5ii MeV 
order=4 'treated as a multiple scatter photon 
go to 21 

end if 


c The charged particles are transported. First the stored parameters 
c are recalled. 


33 


do 120 while (charge .ge. 1) 


e=kinrg (charge) ‘kinetic energy 


cal=kincal(charge) 'particle’s phantom angles 
sal=kinsal (charge) 

cbeta=kincbeta(charge) 

sbeta=kinsbeta(charge) 


xs=kinx (charge) 'particle’s position 
ys=kiny(charge) 
zs=kinz(charge) 


order=kinorder(charge) ‘order of scattering 
nrgint=24 

if (order .eq. 1)then 

chargpart=chargpart+t+l ‘number of charged particles 
end if 


The location and distance to the element wall in the particles direction 
is determined. 


Do only 
depends 


call where 
ro=floaty(phantom(i, jy, k)+1000)/1000. 0 


while the particle’s energy is above the cut-off. The cut-off 
on the density of the phantom element. 


do 110 while (e .gt. O. 1i#ro) 


call where 


c The following statements find the distance to the voxel wall 
c in the coordinate directions. The integer voxel locators are also found. 


nearx=abs(xs—-floaty(i)#diml) 
nearinti=itl 
if (nearx .gt. O. 5#diml) then 
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nearx=diml—-nearx 
nearinti=i-1 
end if 


neary=abs(ys-floaty(j)#diml) 
nearinty=j+1 
if (meary .gt. O. S#diml) then 
neary=diml-neary 
nearint y=j-1 
end if 


nearz=abs(zs—-floaty(k)#dimlz) 

nearintk=k+1 

if (mearz .gt. 0. S¥#dimlz) then 
nearz=dimlz—-nearz 
nearintk=k-1 

end if 


c The nearest voxel is found. The shortest distance to a voxel wall 
c in the coordinate directions, ‘near’, is found 


near=diml 


if (nearx .1t. near) then 
near=nearx 
neari=nearinti 
near j=yJ 
neark=k 

end if 


if (neary .1t. near) then 
near=neary 
neari=i 
nearT jy=nearinty 
neark=k 
end if 


if (nearz .1t. near) then 
near=nearz 
neari=i 
near yj=J 
neark=nearintk 
end if 


c The stopping power, ’s’ is found using log-linear interpolation. 
do 55 while (e .1t. nrg(nrgint)) 
he) nrgint=nrgint-1l 
fracil=(e-nrg(nrgint))/(nrg(nrgint+1)—nrg(nrgint) ) 
s=stop(nrgint)+(stop(nrgint+1)—-stop(nrgint) )#fraci 
s=exp(s) 
ro=floaty(phantom(i, y, k)+1000)/1000. 0 
maxpl=e/s/ro !maximum path length 
if (maxpl .1t. near) go to 112 ‘electron can not leave the voxel 


c The scattering power, ‘th’ is found 


th=tha(nrgint)+(tha(nrgint+1)—-tha(nrgint) )#fraci 
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th=exp(th) 


‘Ep’ is the amount of energy lost in the step. ‘Pl’ is the pathlength 
of the step. 


ep=(1-cay)*e 
pl=ep/s/ro 


Tf ‘pl’ is greater than short, then the pathlength must be shortened 
to ‘short’ so that the particle does not leave the voxel. If so, the 
new energy lost, ‘ep’, must be calculated 


if (pl .gt. short .and. short .gt. 0.0) pl=short 
ep=pl#s#ro 


x=pl#sal¥cbeta 'distance travelled during the step 
y=pl#sal#sbeta 
z=pl#cal 


XxS=xst+x ‘position of the particle at the end of step 
ys=ysty 
zS=zS+z 


oldi=i 'the voxel location before the step is stored 
oldj=y 
oldk=k 
i=yifix(xs/diml)+1 ‘voxel location after the step 
J=yifix(ys/diml)+1 
k=yifix(zs/dimlz)+1 
e=e-ep ‘energy after step 
g=sqrt(thepl#ro)#gaus( yifix(100. O#ran(ili))+1) ‘zenith angle 
r=6. 2832#ran(il) ‘azimuth angle 
Transform from the particle to the phantom coordinate system. 
call geom(r,g,cal, sal, cbheta, sbeta) 
The following statement checks for the edge of the phantom. 
if (xs .le. 0.0001 .or. xs .ge. 20. O#diml) go to 115 
if (ys .le. 0.0001 .or. ys .ge. 20.0#diml) go to 115 
ifieaCzsroles 10.0001 Sand, cal. it. 0:0) goto 115 
if (zs .ge. 20.0#dimlz) go to 115 


By finding the smallest of the old or new voxel location the correct 
voxel to assign the dose is guaranteedi even if the particle is on 
a voxel wall. 


smalli=jminO(oldi, i) 
smally=jminO(oldy, Jj) 
smallk=yjyminO(oldk, k) 


Assign the energy lost during the step as dose. 
ro=floaty(phantom(smalli, smally, smallk)+1000)/1000. 0 


dose(order, smalli, smally, smallk)= 
1 dose(order, smalli, smally, smallk)+ep/ro 
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110 continue 'charged particle transport step ends here 


Cc 
Cc 


n 


Cc 
c 
c 


The following checks to see if the nearest voxel is outside the phantom. 
Tf it is half of the remaining energy will be deposited as dose 


if (meari .le. O .or. neari .ge. 21) then 
e=e/2.0 
go to il2 

end if 


if (neary .le. O .or. neary .ge. 21) then 
e=e/2.0 
go to i112 

end if 


if (meark .le. O .or. neark .ge. 21) then 
e=e/2.0 
go to i112 

end if 


‘Nearo’ is the density of the nearest voxel. The remaining energy is 
assigned to the last and nearest voxel such that they both get equal 
energies. 


nearo=floaty(phantom(neari, neary, neark)+1000)/1000. 0 
dose(order, i, jy, k)=dose(order, i» jy, k)+e/(rotnearo) 
dose(order, neari, neary, neark)= 

1 dose(order, neari, near y, neark)+e/(ro+nearo) 
ga to 113 


Statement 112 and its following statement are only done if it is found 
that the particle can not leave the phantom or if the neighboring voxel 
is outside the voxel. 


112 ro=floaty(phantom(i. jy, k)+1000)/1000. 0 


dose(order, i, y,k)=dose(order, i. jy, k)+e/ro 


py Re continue 


c 


The following statements are done if the charged particle was a positron. 
if (kintype(charge) .eq. -1) then 


kintype(charge)=0 
charge=charge-1 
annhilflag=1 'sets flag for the second photon 


ri=ran(il) 

cal=cos(6. 2832#ri1) 
sal=sqrt(1. O-cal#cal) 
r2=ran(il) 
cheta=cos(6. 2832#r2) 
sbeta=sin(6. 2832*r2) 


nrgx=0. Sil 
photox=xs ‘the last position of the 
photoy=ys ‘positron must 


photoz=zs 'be stored 
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order=4 ‘include annhilation photons as 
‘multiple scatter 


go to 21 
end if 
LS continue 


charge=charge-1 


120 continue tend of loop for the charged particle histories 
180 continue ‘end of loop for the photon histories 
200 continue ‘end of loop for the spectral bins 


c The following statements normalize the dose to the maximum 
c dose encountered and writes the normalized dose 


total=0.0 
do 205 k=1,20 
do 205 j=1, 20 
do 205 i=1,20 
205 total=total+dose(i,i, j,k) ‘total primary KERMA 


c Print the beam parameters. 


write(4,#) ‘Spectrum Multiplier=’, multi 

write(4,#) ‘# Photon Histories=’, history 

write(4,%*#) ‘Shiftx=’, shiftx 

write(4,#) ‘’Shifty=’, shifty 

write(4,#) ‘Primary Interaction Depth Integer=’, intest 
write(4,%) ‘X-Y Element Size=’, diml 

write(4,*) ‘Depth Dimension=’, dimlz 

write(4,#) ‘Total Primary Dose=’, total 

write(4,*) ‘Number Of Primary Charged Particles’, chargpart 


do 210 k=1,25 
do 210 j=1,25 
do 210 i=1,25 
do 210 q=1i,5 
doset(i, jy, k)=doset(i, jy, k)+dose(q, i, yj, k) 'total dose 
‘distribution 


210 continue 


c Two-dimensional primary dose distribution in the plane of the 
c primary pencil beam. The dose is normalized to the total primary KERMA. 


write(4,#) ‘Coronal Primary Dose Slice’ 
do 215 k=1,20 
do 212 i=1,20 
212 atr(i)=dose(1,i,11,k)/total 
215 write (4,650) (arr(i),i=5, 17) 


c Three-dimensional total dose distribution normalized to the primary KERMA. 


do 229 k=1,20 
write(4,%#)’ “ 
write(4,580) k 


aan i 
rortelidane siubons* ; : 
Ta sShze: WL QEP hm i> 
2 = ; 
yeh oe y ni dnes ee 
oe A : snes 
Ne Ste | 
/ : 
hearai> eft) Tot Quel te love 
nosvoiq =af 40% deat 4a? geves 
4 | 


a y22dn4, o29 Fo? Rese Va, tome 


2 ) | a : Bit 
momirem, ot? ooeeae ug ties sama % ett ime 
each badis aan tite ore ; oe a: . 


4 ‘ ? ; = a 
to anaiq a fr! isi Lae vbr 
igjot aft ps wens isd ain aonb 
exile ueng i 
eye 
Nees abe 
oni < oe 


425 


do 225 j=1, 20 
do 224 i=1,20 
224 atr(id=doset(i, jy, k)/total 
225 write(4,650) (arr(i), i=5,17) 


c Three-dimensional distributions of the primary KERMA 
c normalized dose components. 


do 229 q=1,4 
write(4,#) ’ ” 
if (q .eq. 1)then 
write(4,.%) ‘Primary Dose’ 
else if (q .eq. 2) then 
write(4,#) ‘First Scatter Dose’ 
else if (q .eq. 3) then 
write(4,#) ‘Second Scatter Dose’ 
else if (q .eq. 4) then 
write(4,#) ‘Multiple Scatter Dose’ 
end if 
do 228 j=1,20 
do 227 i=1,20 


227 atr(id=dose(q, i, jy, k)/total 
228 write(4,650) (arr(i), i=5, 17) 
229 continue 


c Three-dimensional distribution of the primary KERMA 
c normalized total scatter dose 


write(4,%) ‘Total Scatter Dose’ 
do 245 k=1, 20 
write (4,%)’ ¢% 
write(4,580) k 
do 245 j=1,20 
do 244 i=1,20 
244 arr(id=(doset(i, jy» k)-dose(l, i, y, k))/total 
245 write(4,650) (arr(i), i=5,17) 


c Three-dimensional distribution of the primary KERMA 
c normalized primary plus first scatter dose 


write(4,#) ‘Primary + First’ 
do 247 k=1,20 
write (4,%)’ % 
write(4,580) k 
do 247 j=1,20 
do 246 i=1,20 


246 artr(i)=(dose(1, i, y, k)+dose(2, i, jy» k))/total 

247 write(4,650) (arr(i), i=5, 17) 

250 format(’ ’,13,2x, #6. 2) 

340 format(’ ’,f5.2,2x,f5. 2) 

345 format(’ ‘’,i2) 

350 format(’ ’,10(f4.2,2x)) 

420 format(’ ’,e10. 4) 

422 format(eld. 4) 

550 format(’ %,10(1x, #6. 4)) 

970 format(’ % €5.2,1x,3¢(F6. 4,1x),2¢F6. 3, 1x)) 

580 FOTMAt ( / HHH EHH EHH HHH KH, 1D, OHHH HHH HHH / ) 
650 format(’ ’,20(1x,e8.3)) 

655 format(’ ‘,1x,al,3x, a4, 4x, a14, 2x, a7, 4x, a5, 3x, a6, 2x, aB) 


7350 format(’ ‘’,1i12,2x, 5. 3,8x, #4. 3, 6x,4(4x, #5. 3) ) 
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close(1) 
close(4) 
close(é6) 
close(8) 
close(ili) 
close(12) 


stop 
end 


subroutine geom(r,g,cal, sal, cheta, sbeta) 


This subroutine performs the rotational transformations between 
the coordinate system ‘attached’ to the particle and the system 
‘attached’ to the phantom. 


real r,g,cal, sal, cbheta, sbeta 
real ss,cc,oldcal, oldsal, oldcbeta, oldsbeta 


‘G’ and ‘r’ are the sine and cosine respectively in the particle system. 
‘Sal’ and ‘cal’ are the sine and cosine respectively of the zenith 
angle in the phantom system. ‘’Sbeta’ and cbheta are the sine and 

cosine respectively of the azimuth angle in the phantom system. 


if (g .eq. 0.0) go to 400 


oldsal=sal 
oldcal=cal 
oldcbeta=cbeta 
oldsbeta=sbeta 


cal=cal*cos(g)-sal¥sin(g)#cos(r) 


if (abs(cal) .ge. 1.0) then 
cal=1. O#sign(1l.0,cal) 
sal=0.0 
cbheta=cos(r) 
sbeta=sin(r) 
go to 400 

end if 


sal=sqrt(i. O-cal#cal) 


if (oldsal .eq. 0.0) then 
cheta=cos(r) 
sbeta=sin(r) 
go to 400 

end if 


ss=sin(g)#sin(r)/sal 
cc=(sin(g)#cos(r)#oldcal+cos(g)#oldsal)/sal 


if (abs(ss) .ge. 1.0) then 
ss=1. O#sign(1.0,s5) 
end if 


if (abs(cc) .ge. 1.0) then 
cc=1. O¥sign(1.0, cc) 
end if 
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cheta=cc#oldcbeta-sstoldsbeta 
sbeta=ss*toldcbetat+cc#oldsbeta 


if (abs(cbeta) .ge. 1.0) cheta=1. O#sign(i. 0, cheta) 
if (abs(sbeta) .ge. 1.0) sbeta=1. O#sign(1. O. sbeta) 


continue 


return 
end 


subroutine where 


This subroutine finds the position and distance to the voxel boundary 


in the particle direction. 


real xs,ys,zs, short, diml, dimlz 
real cbheta, sbeta, cal 

real lx,ly,1lz 

integer i, yk, ib» yb, kb 


common xs,ys, 2s, i1, jy, k, cbheta, sbeta, cal, sal. diml, dimlz, short 


Determine 


Determine 
direction 


which element the photon is in. 


i=yifix(xs/diml)+1 
J=sifix(ys/diml)+1 
k=yifix(zs/dimlz)+1 


the position of the voxel’s walls in the 


in which the pazrticle is headed. 


if (sal .eq. 0.0 .or. cbeta .eq. 0.0) 
1x=999. 9#diml 
go to 10 

end if 


if (cbeta .1t. 0.0) then 
if (xs .eq. (floaty(i-1))#diml) then 
lx=(floaty(i-2)#diml-xs)/sal/cbheta 
else 
lx=(floaty(i-1)#diml-xs)/sal/cbeta 
end if 
else 
lx=(floaty(id#diml-xs)/sal/cbheta 
end if 


if (sal .eq. 0.0 .or. sbeta .eq. 0.0) 
1y=999. 9#diml 
go to 20 

end if 


if (sbeta .1t. 0.0) then 
if (ys .eq. (floaty(y-1))#diml) then 
ly=(floaty(y-2)*#diml-ys)/sal/sbeta 
else 
ly=(floaty(j-1)*diml-ys)/sal/sbeta 
end if 
else 
ly=(floaty(j)#diml—-ys)/sal/sbeta 
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end if 
20 if (cal .eq. 0.0) then 
12=999. G#edimlz 
go to 30 
end if 


if Ccal sit. O20) then 
if (zs .eq. (floaty(k-1))#dim1lz) then 
lz=(floaty(k-2)#dimlz—-zs)/cal 
else 
lz=(floaty(k-1)#dimlz-zs)/cal 
end if 
else 
lz=(floaty(k)#dimlz—-zs)/cal 
end if 


c Determine the distance to the wall that the particle 
c in the direction of the particle 


30 short=1z 
if (ly .1t. short) short=ly 
if (lx .1t. short) short=l1x 


return 
end 
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Appendix ii. Listing of the program Voive. for 


NMnNnANnnNnnnNnnnaAnnan 


Program Volve 


This program convolves the fluence with dose spread arrays to 
yield a dose distribution. It is carried out in "the interaction 
point of view". The primary and truncated first scatter (TFS) 
dose spread arrays are combined. This fine resolution dose spread 
array is is convolved first. Then a lower resolution, residual 
first and multiple scatter dose spread array is convolved. The - 
homogeneous and heterogeneous dose distribution is calculated 
separately and a correction factor is determined at the end. The 
dose due to contamination is included. The dose distribution 

can be smoothed. The variables are documented where they occur. 
Lines which have been "commented out" are optional. 


integer albnd, aubnd, blbnd., bubnd, depth, lbndi, ubndi 
integer albnd2, aubnda2, blbnd2, bubnd2, depth2 
integer lbndi2, ubndi2 

integer lbndy, ubndy, lbndk, ubndk 

integer lbndy2, ubndy2, lbndka, vbndk2 

integer widthi, heighti, lengthi 

integer i, j,k, ii, yy, kk» deli, dely, delk 

integer clipwid,cliphei,cliplen, roe, roe2, kprime 


real mu, ro(-25: 25, -25: 25, 50), ro2(-5: 5, -5: 5, 10) 

real invsqr, ssd,dref 

real spray1(5, —-6: 6, -6: 6, -2: 12) 

real rfmspray(5,; -6: 6, -6: 6, -2: 12) 

real fluxinc, fluxfact 

real avdens, avdensd, avspray 

real dosehet(-25: 25, -25: 25, 50), dosehet2(-5: 5,-5:5,0:10) 
real dosehom(-25: 25, —-25: 25, 50), dosehom2(-5: 5, -5:5,0:10) 
real kerma(50), kerma2(10) 

real oldkerma, cf(-25: 25, -25: 25, 50), kermhet, kermhom 
real arrcf(20), arrhom(20), arrhet(20) 

real £0, fl, f2,c¢c1,¢c2, row 

real hom(-25: 25, -25: 25, 50), het (-25: 25, -25: 25, 50) 

real contam(0: 50,2), avcontam(2), resolv 


common deli,dely,delk, i, jy, k, ro, avdens 


Spray. dat contains the primary plus truncated first scatter (TFS) 
dose spread arrays. 


Cf. dat contains the output homogeneous and heterogeneous dose 
distributions and the correction factor distribution. 


Rfmspray.dat contains the residual first and multiple scatter (RFMS) 
dose spread arrays. 


open(unit=1, status=’old’, file=’spray. dat’) 
open(unit=3, status=‘’new’, file=’cf. dat’) 
open(unit=4, status=’old’, file=’rfmspray. dat’) 


User requested statements. A slab geometry is modelled by a cubic 
water-like phantom. The beam is parallel. 


ssd=100. 0 'source to surface distance 
dref=3. 5 'dmax 
resolv=1.0 ‘specifies voxel size for contamination calculation 
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2-Dimensional 


albnd=-4 'beam half-width in negative x-direction 
aubnd=4 'heam half-width in positve x-direction 
bibnd=-4 ‘beam half-width in negative y-direction 
bubnd=4 ‘beam half-width in positive y-direction 
depth=20 ‘phantom depth 


3-Dimensional 


The calculation windows may be 1,2 or 3-dimensional. 


lbndi=albnd ‘calculation window in negative x-direction 
ubndi=aubnd ‘calculation window in positive x-direction 
lbnd y=blbnd ‘calculation window in negative y-direction 
ubnd y=bubnd ‘calculation window in positive y-direction 


lbndk=1 
ubndk=depth 


lbndi=albnd-1 
ubndi=aubnd+1 


‘lower calculation window in z-direction 
‘upper calculation window in z-direction 


(in the plane of the 


‘calculation 
‘calculation 


central axis) 


window 
window 


in 
in 


negative 
pasitive 


x-direction 
x-direction 


lbnd y=0 ‘calculation window in negative y-direction 
ubnd j=0 ‘calculation window in positive y-direction 
lbndk=1 ‘lower calculation window in z-direction 


ubndk=depth 


‘upper calculation window in z-direction 
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1-Dimensional (off axis at a specified depth) 


lbndi=albnd-1 
ubndi=aubnd+1 


‘calculation window in negative x-direction 


NmnNnnaAnn 


n 


1-Dimensional 


‘calculation window in positive x-direction 


lbnd jy=0 ‘calculation window in negative y-direction 
ubnd j=0 ‘calculation window in positive y-direction 
lbndk=16 ‘lower calculation window in z-direction 
ubndk=16 ‘upper calculation window in z-direction 


(along the central axis) 


lbndi=0 ‘calculation window in negative x-direction 
ubndi=0 ‘calculation window in positive x-direction 
lbnd y=0 ‘calculation window in negative y-direction 
ubnd jy=0 ‘calculation window in positive y-direction 
lbndk=1 ‘lower calculation window in z-direction 


ubndk=depth 


‘upper calculation window in z-direction 


clipwid=2 'a parameter for delimiting width of dose spread array 
cliphei=0 'delimits height (delk=negative) 
cliplen=0 'delimits length (delk=positive) 


Read statements. Read height, length, width, and values of the 


dose spread arrays, and the density array. 


read(1,3000) heightil, lengthi, widthil 
do 100 n=1,5 
read(1, #) 
read(4, #) 
do 100 k=—-height1, lengthi 
read(1,*) 
read(1, *) 
read (4, #) 
read (4, *) 
do 100 jy=-widthi,widthi 
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read(1,3010) (sprayi(n, i, jy. kk), i=-widthi, widthi) 
read(4,3010) (rfmspray(n, i, jy. k), i=-widthi, widthi) 
continue 


read (4, #) 
read (4, #) 
do 125 n=1,5 

read (4, *) 

read (4, #) 

do 125 k=-1,3 

read (4, *) 

do 125 j=-1,1 

Tead (4,3010) (rfmspray(n, i. yk), i=—-1, 1) 

continue 


c The relative fluence in a homogeneous 
c phantom as a function of depth is calculated 


150 


Nnann 


kerma(1)=1. 00 

mu=0. O31 

atten=exp (-mu) 

do 150 k=2, depth+1 
kerma(k)=kerma(k-1)#atten 

continue 


"Contam(depth, region)" is the measured depth dependence of contamination. 
The depth is measured in centimeters. Region=1 is inside the field 
Region=2 is outside the field. See Volume 1 and Equation 7.4.9 in 

Volume 2. 


contam(O, 1)=1. 00 
contam(1i,1)=0. 89 
contam(2, 1)=0. 50 
contam(3, 1)=0. 28 
contam(4, 1)=0. 18 
contam(5, 1)=0. 12 
contam(é6, 1)=0. 10 


contam(0O, 2)=1. 00 
contam(1,2)=0. 94 
contam(2, 2)=0. 64 
contam(3, 2)=0. 44 
contam(4, 2)=0. 37 
contam(5, 2)=0. 33 
contam(é6, 2)=0. 30 


c Beyond dmax the contamination decreases exponentially. 


175 


atten=exp(-0. 050) 
do 175 n=1,.2 
do 175 k=7,50 
contam(k,n)=contam(k—-1,n)#atten 
continue 


c The following statements describe a slab phantom. The central slab can be 
c made up of two regions 


do 200 k=1,15 
do 200 j=-25,25 
do 200 i=-25, 25 
To(i, je kK)=1.0 
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200 continue 


do 220 k=16,23 
do 220 jy=-25,25 
do 220 i=-25,.3 
ro(is j,kd)=0.3 
220 continue 


do 230 k=16,23 
do 230 j=-25,25 
do 230 i=4, 25 
ro (i, yj» k)=0.3 
230 continue 


do 240 k=24, 50 
do 240 j=-25, 25 
do 240 i=-25, 25 
Toi, Je k)=1.0 
240 continue 


c Convolve for all pencil beams at all depths by summing over all primary 
c interaction voxels. "Fluxfact" is the ratio of heterogeneous to homo-— 
c geneous fluence at a point in a pencil beam. 


do 600 i=albnd-1, aubnd+1 ‘increment in the x-direction 
ii=ynint(floaty(i)/5. 0) 'for RFMS in the x-direction 
do 600 j=blbnd-1, bubnd+1 ‘increment in the y-direction 
Js=ygnint(floaty(y)/5. 0) 'for RFMS in the y-direction 


c The pencil beams just outside the field can be given a smaller fluence 
c This can account for a geometrical penumbra. ‘Fluxinc is the incident 
c relative fluence. 


Cc Wf (i eq. albnd—1 jor.) 2) eq. aubnd+!l or: 

c 1 j eq. bi bnd=tenor yo eq.  bubndti ) then 
c fluxinc=0. 20 

Cc else 

c fluxinc=1. 00 

c end if 


c The following can describes a bar shield. 


G if (ie eqn 1) Ores le eq On onre ee tereds) ao) then 
c fluxinc=0. 031 
Cc end if 
fluxfact=1. 00 'fluxfact is the fluence compared to a 
‘homogeneous phantom 
do 600 k=1, depth ‘increment depth (z-direction) 
kKk=(k-1)/5+1 'for RFMS in the z-direction 


c An inverse square factor can reduce the primary fluence. 


c invsqr=(ssdt+dref) ##2/ 
(2 xf ((ssdt+floaty(k)—-. 5) ##2+ 
c 2 (floaty (yj) ##2+( floaty (i) )##2) 


invsqr=1.0 
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"Kermhom" and "“kermhet" are factors dependent only on the point 
of interaction of the primary fluence 


kermhom=kerma(k)#invsqr#fluxince 
kermhet=kermhom#ro(i, y, k)#fluxfact 


Sum over all contributions of the fine resolution dose spread array 
by summing over all dose deposition voxels 


do 500 deli=-widthit+clipwid., widthi-clipwid 
do 500 dely=-widthi+tclipwid, widthi-clipwid 
do 500 delk=—-heighti+cliphei., lengthi-cliplen 


If the contribution is outside the dose calculation window 
don’t bother summing. 


Lee Citdel2 elt lel binds). gio. to .500 
if (itdeli .gt. ubndi) go to 500 
if Cytdel) lt. lbndy) go to S00 
Lt Cy tdet meng tamuUD nia) goto .cOO 
if (k+delk .1t. lbndk) go to 500 
if (kt+tdelk .gt. ubndk) go to 500 


Subroutine "averdens" returns "“avdens", the average density 
between the primary interaction and dose deposition voxels. 
"Roe" is an integer identifying a dose spread array obtained for 
a specific density (roe=1 for a physical density of 0.2 g/cm##3, 
Toe=2 for 0.4 g/cm##3 ... roe=5 for 1.0 g/cm##3) 


call averdens 


"Avspray" is the dose/KERMA value linearly interpolated between dose 


spread arrays at fixed densities 


avspray=sprayl(roe, deli,dely,delk)+ 
(sprayl(roe+i, deli, dely, delk)—- 
sprayl(roe, deli,dely, delk))+* 
(5. Otavdens—floaty (roe) >) 


W Me 


Second order interpolation can be used if desired. 


if (avdens .eq. 1.0) then 
avspray=sprayi(5,deli,dely, delk) 
else 
if (avdens .ge. 0.6) then 
fO=sprayl(roe-1i, deli, dely, delk) 
fl=sprayi(roe, deli, dely, delk) 
f2=spraylil(roeti, deli, dely,delk) 
Tow=(avdens*5. O-floaty(roe))/5. 0+0. 2 
else 
#O0=sprayi(roe, deli, dely, delk) 
Fi=sprayi(roe+i, deli, dely, delk) 
f2=sprayll(roeta2, deli, dely, delk) 
row=(avdens*#5. O-floaty(roe))/5.0 
end if 


ci=(row/0. 2)#(-3. OF FO+4. OX F1-F2)/2. 0 
c2=(row/0. 2) ##2#( FO-2. OX F1l+F2)/2. 0 


avspray=f0+cl+c2 
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c end if 
c "Dosehet" is the heterogeneous dose 
dosehet(itdeli, y+dely, k+delk)= 
1 dosehet(it+deli, j+dely, k+delk)+ 
2 kermhet/avdens#avspray 
c "Dosehom" is the homogeneous dose 
dosehom(i+deli, ytdely, ktdelk)= 
1 dosehom(i+deli, y+dely, kt+delk)+ 
2 kermhom#spray1(5,deli,dely,delk) 


500 continue 


fluxfact=fluxfact#exp (mu#(1. O-ro¢(i,s y»k))) !cumulative relative fluence 


avdens2=avdensa2+ro(i, jy, k) 'density sum for RFMS 
Toglii, yy kk=roQlii, yy, kkd+roci.s yr k?) 'for RFMS average voxel density 
600 continue 


c At this point the primary plus truncated first scatter dose spread arrays 
c have been convolved. The residual first and multiple scatter (RFMS) dose 
c spread array is now convolved. 


albnd2=yjnint(floaty(albnd-1)/5. 0) 'field boundaries 
aubnd2=ynint(floaty(aubnd+1)/5. 0) 
blbnd2=jnint(floaty(blbnd—-1)/5. 0) 
bubnd2=jnint(floaty(bubnd+1)/5. 0) 

depth2=(depth-1)/5. O+1 


lbndid=ynint(floaty(lbndi)/5. 0)-1 'dose calculation window 
ubndid=jnint(floaty(ubndi)/5. O)+1 
lbndy2=snint(floaty(lbndj)/5.0)-1 
ubndy2=ynint(floaty(ubndy)/5. O)+1 

lbndk2=0 

ubndk2=(ubndk-1)/5+2 


c Calculation of fluence ina homogeneous phantom taking into account 
c the coarser resolution. 


kerma2(1)=1. 00 
do 700 k=2, depth2+1 
atten=exp(—mu*5. 0) 
kerma2(k)=kerma2(k—-1)#atten 
700 continue 
c Calculate the average density in the irradiated portion of the phantom. 


avdens2=avdens2/ 
1 ((aubnd-albnd+3)#(bubnd—-blbnd+3)#depth ) 


roe=jifix(avdens2#s5. 0) 
roe2=(5. O#xavdens2-floaty(roe) ) 


do 1200 i=albnd2, aubnd2 ‘increment in x-direction 
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do 1200 y=blbnd2, bubnd2 ‘increment in y-direction 


fluxinc=1. 00 


n 


Calculate the fluence near the field boundary if the geometrical 
penumbra is being taken into account. 


n 


c if (i .eq. aubnd2) then 

c fluxinc=(floaty(jyiabs(aubnd)+3) 

c 1 73. O-floaty(jyiabs(i)d))#fluxince 
c else if (i .eq. albnd2) then 

c fluxinc=(floaty(jyiabs(albnd)+3) 

(e 1 735. O-floaty( jiabs(i)))#fluxinc 
(= end if 

c 

c if (jy .eq. bubnd2) then 

( fluxinc=(floaty( jiabs(bubnd)+3) 

(= 1 73. O-floaty(jiabs(j)))#fluxine 
c else if (jy .eq. blbnd2) then 

Cc Ffluxinc=(floaty(jyiabs(blbnd)+3) 

Cc 1 7/3. O-floaty(jiabs(y)))#fluxine 
c end if 

c The following is used with the bar shield 

(= if (i .eq. O) then 

c fluxinc=0. 40 

c end if 


c The following is used with a wedge 


c fluxinc=fluxinc#O. 5#floaty(i-albnd2)/ 
c 1 float y(aubnd2-albnd2)+0. 5#fluxinc 
do 1200 k=1,depth2 ‘increment depth (z-direction) along 


‘pencil beam 


c Used with the inclusion of inverse square attenuation. 


c invsqr=(ssd+dref) ##2/ 
c 1 ((ssd+5. O#floaty(k)-2. 5) ##2+ 
c 2 (5. O#Floaty(y) ##2+(5. O#Floaty(i) )##2) 


invsqr=1. 00 


roa(i, y, k)=roa(i, jy, k)/125. 0 'finish average voxel density 
kermhom=kerma2(k)#invsqr*fluxinc ‘homogeneous fluence 
kermhet=kermhom#ro2(i, jy, k) ‘heterogeneous fluence*voxel density 


c Sum over all contributions of the RFMS dose spread array by summing 
c over all dose deposition voxels. 


do 1100 deli=-widthi, widthi 
do 1100 dely=-widthi,widthi 
do 1100 delk=—-heighti1, lengthi 


if (it+deli .1t. lbndi2) go to 1100 ‘inside dose window? 
if (itdeli .gt. ubndi2) go to 1100 
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c Interpolate 


1100 
1200 


ce Print 


1201 


1202 


1203 


1204 


ifeuC jeden ge cet. 
if (ytdely .gt. 
if (k+delk .1t. 
if (k+delk . gt. 


avspray=rfmspray (roe, deli,dely, delk)+ 
1 (rfmspray(roet+i, deli, del yj, delk)-— 


lbnd j2) 
ubnd ya) 
lbndk2) 
ubndk2) 


go 
go 
go 
go 


to 
to 
to 
to 


the heterogeneous dose spread 


1100 
1100 
1100 
1100 


array. 
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2 rfmspray(roe, deli,dely, delk))#roe2 


dosehet2(itdeli, ytdely, ktdelk)= 
1 dosehet2(itdeli., y+dely, k+delk)+ 
2 kermhet/avdens2tavspray 


dosehom2(i+deli, ytdely, k+delk)= 


at dosehom2a(i+deli, y+dely, kt+delk)+ 

2 kermhom#rfmspray(5,deli,dely, delk) 
continue 

continue 


the primary plus TFS dose before smoothing. 


write(3,+#) “Homogeneous Dose Before Smoothing 


do 1202 k=lbndk,. ubndk 


write(S, #) “wee eH KH, ky, RHEE 


do 1201 y=lbndy,ubndy 
do 1201 i=lbndi,ubnd 

atrhom(i)=dosehom(i 
continue 


1 
» Je k) 


write(3, 3100) (arrhom(i), i=lbndi, ubndi) 


continue 


‘increment hetero dose 


‘increment homo dose 


(No Multiple) ’ 


write(3, *#) ‘Heterogeneous Dose Before Smoothing (No Multiple)’ 


do 1204 k=lbndk, ubndk 


write (3, #) “Hee eHeHHH KH, kk, “HHH HHH 


do 1203 j=lbndjy, ubndy 
do 1203 i=lbndi, ubnd 
arrhet(i)=dosehet(i 
continue 


1 
rade Kad 


write(3, 3100) (arrhet(i), i=lbndi, ubndi) 


continue 


c "9-point" smoothing. 
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delyj=0 

do 1210 k=lbndk+1,ubnd 
do 1210 i=lbndi+i1, ub 


k-1 
ndi-1 


do 1210 deli=i-1, iti 


do 1210 delk=k-1, 


k+1 


het(i, yp» k)=het (i, y, k)+dosehet(deli,dely, delk) 


continue 


do 1220 k=lbndk+1, ubndk—-1 
do 1220 i=lbndi+t+1, ubndi-1 


dosehet(i, jy, k)=het(i, yo kd/9.0 
dosehom(i, y, k)=hom(is y,k)/9.0 


hom(i, y, k)=hom(i., y, k)+dosehom(deli, dely, delk) 
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c het(i, Je kd» =0.0 
(Ss hom(i, Je k)=0.0 
1220 continue 


c Print the primary plus TFS dose after smoothing 


c write(3,#) ‘Homogeneous Dose After Smoothing (No Multiple) ‘ 
c do 1250 k=lbndk, ubndk 

Cc write (3, +) “wee k=, kK, (eH 

c do 1250 y=lbndy,.ubndy 

(= do 1225 i=lbndi, ubndi 

c arrhom(i)=dosehom(i, y, k) 

1225 continue 

( write (3, 3100) (arrhom(i?), i=lbndi, ubndi) 

1250 continue 


c write(3,#) ‘Heterogeneous Dose After Smoothing (No Multiple)’ 
c do 1300 k=lbndk, ubndk 

c write (3, #) (“eee kK’, k, HHH HHH / 

( do 1300 j=lbndj, ubndy 

c do 1275 i=lbndi, ubndi 

c arrhet(i)=dosehet(i, jy, k) 

1275 continue 

c write (3, 3100) (arrhet(i), i=lbndi, ubndi) 

1300 continue 


c Print the phantom and voxel densities used for the RFMS convolution. 
Cc write(3, *) ‘Average Phantom Density’, avdens2 


c write(3, *) ‘Average Voxel Density’ 

c do 1310 k=lbndk2+1, ubndk2-1 

c write (3, *) “HeHHHHHH KE, kK, HHH HHH / 
Cc do 1310 y=lbnd y2+i, uvbnd y2-1 

c do 1305 i=lbndi2+1, ubndi2-1i 

(= arrhom(i)=ro2 (i, jy, k?) 


£305) continue 
¢ write(3, 3100) (arrhom(i), i=lbndi2+1, ubndi2-i) 
1310 continue 


c Print the RFMS dose distribution. 


c write(3,#) ‘Homogeneous Multiple Dose’ 
Cc do 1350 k=lbndk2, ubndk2 

c write (GS, #) ‘eee k=! kk, HHH HHH 
Cc do 1350 jy=lbndy2, ubndy2 

c do 1325 i=lbndi2, ubndi2 

Cc atrhom(i)=dosehoma(i, jy, k) 


1325 continue 

c write(3, 3100) (arrhom(i), i=lbndi2, ubndi2) 
1350 continue 

c write(3, #) ‘Heterogeneous Multiple Dose’ 

c do 1400 k=lbndk2, ubndk2 

(2 write (GS, #) “Hee HH KH, k, “HHH “ 

(e do 1400 j=lbndj2, ubndy2 

c do 1375 i=lbndi2, ubndi2 

Cc arrhet(i)d=dosehet2(i» y, k) 

1375 continue 

c write (3, 3100) (arrhet(i)d, i=lbndi2, ubndi2) 


1400 continue 
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c Add the RFMS dose distribution on to the primary plus RFS dose 
c distribution. Calculate the radiological depth for inclusion 
c of the contamination dose. 


1500 


1550 


j=0 
del y=0 
do 1500 k=lbndk, ubndk 
do 1500 i=l1bndi, ubndi 
do 1500 deli=i-2. i+2 
do 1500 delk=k-—2, k+2 
het (i, jy, k)=het(i, yk) 
+dosehet2( ynint(floaty(deli)/5.0) 
»jnint(floaty(dely)/5. 0) 
» (delk-1)/5+1) 


QMUre 


hom(i. ys k)=hom(i, y, k) 
1 +dosehom2( ynint(floaty(deli)/5. 0) 
2 1 Jnint(floaty(dely)/5. 0) 
Ss » (delk-1)/5+1) 


continue 
do 1575 i=lbndi, ubndi 
avdens=0. 0 


do 1575 k=lbndk, ubndk 


dosehet(i, jy, k)=het(i, jy, k)/25. O+dosehet (i, j,k) 
dosehom(i, y, k)=hom(i, jy, k)/25. O+dasehom(is y, k) 


avdens=avdens+ro(i, j,k?) ‘radiological depth 


do 1550 n=1i,2 


avcontam(n)=contam( jyifix(avdens),n)+ ‘calculate the 
1 (contam( yifix(avdens)+i,n)- 'contamintion 
2 contam(jifix(avdens),n))+* 'at the radiological 
3 (avdens-floaty(yifix(avdens))) 'depth 


continue 


c Calculate the contamination dose. The contamination dose has a Gaussian 


c dependence on position in the field from the central axis, a linear increase 


c with field size and a depth dependence found from a look-up table 
if (i .ge. albnd .or. i .le. aubnd) then 


dosehet(i, y, k)=dosehet(i, y, k)+ 


1 avcontam(1)#sqrt(floaty((aubnd-albnd)+# 
2 (bubnd—-blbnd)))#0. Ol#resolv 
dosehom(i, jy, k)=dosehom(i, y, k)+ 
1 contam(k, 1)#*sqrt(floaty((aubnd-albnd)+# 
2 (bubnd—blbnd)))+#0. O1 
else 


dosehet(i, y, k)=dosehet(i, jy, k)+ 
i avcontam(2)#sqrt(floaty((aubnd-albnd)+# 
2 (bubnd—-blbnd)))#0. Ol#resolv# 
3 exp(-2. O#floaty(i*#i)/ (float j(aubnd-albnd) )##2)+# 
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4 exp(—-2. O¥floaty(y#y)/(Ffloaty(bubnd-blbnd) )##2) 


dosehom(i, jy, k)=dosehom(i, y,k)+ 


1 contam(k, 2)*sqrt(floaty((aubnd-albnd)# 

2 (bubnd—-blbnd)))#0. 01% 

3 exp (—-2. O#floaty(i*id/(floaty(aubnd-albnd) )##2)* 

4 exp (-2. O#floaty( yey)/ (floaty (bubnd—blbnd) )##2) 
end if 

continue 


c Calculate and print the correction factor. 


1700 


1800 


1900 


2000 


2100 


2200 


3000 
3010 
3020 
3060 
3070 
3100 


write(3,*) ‘Correction Factor’ 
do 1800 k=lbndk, ubndk 
write (3, *) “Kee KHL kK, (RRR 
do 1800 j=lbndy, ubndy 
do 1700 i=l1bndi, vbndi 
cf(i, y,k)=dosehet(i, jy, k)/dosehom(i, j,k) 
arrcf(id=cf li, yk) 
atrhom(i)=dosehom(i, J. k) 
arrhet(id=dosehet(i, j,k) 
continue 
write(3, 3100) (arrcf(i), i=lbndi, ubndi>) 
continue 


write(3, #)’Total Smoothed Homogeneous Dose’ 
do 2000 k=lbndk, ubndk 
write (3, ) RRR HHH KH, ky HHH” 
do 2000 j=lbndy,ubndy 
do 1900 i=lbndi, ubndi 
arrhom(i)=dosehom(i, jy, k) 
continue 
write(3, 3100) (arrhom(i), i=lbndi, ubndi) 
continue 


write(3,#) ‘Total Smoothed Heterogeneous Dose’ 
do 2200 k=lbndk, ubndk 
write (3, +) “wee HHH KH, kK, HRA HEH 
do 2200 j=lbndy,ubndy 
do 2100 i=lbndi, ubndi 
arrhet(i)d=dosehet(i, j,k) 
continue 
write(3, 3100) (arrhet(i), i=lbndi, ubndi) 
continue 


format(’ ’,3¢12,1x)) 
format(’ ’,25(#8. 6,1x)) 
format (25(#4. 2,1x)) 
format(f5. 1) 

format (’ ’, #46. 3) 
format(’ ’,25(f5.3,1x)) 


stop 
end 


subroutine averdens 


integral differential analyzer algorithm to find the 


average 
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c density along a path between the primary interaction and dose 
c deposition voxels. The algorithm samples the density along 
¢ regular increments along the path 


integer i, y,k-deli,dely,delk, longest 


real ro(—-25: 25, -25: 25, 50), avdens, xinc, yinc, zine 
real scale, denom 


common deli,dely,delk, i, jy. k, ro, avdens 


if (deli .eq. O .and. dely .eq. O .and. delk .eq. O) then 
avdens=ro(i, j,k) 

go to 200 
end if 


c The longest dimension between the interaction and dose deposition 
c sites is found. 


longest=jmaxOlabs(deli), abs(dely), abs(delk)) 
c The increments in each direction between the sample points. 


xinc=floaty(deli)/floaty(longest) 
yinc=floaty(deljy)/floaty(longest) 
zinc=floaty(delk)/floaty(longest) 


c The weight of the interaction voxel in the average density can 
c be varied for different dose deposition voxel locations. 


if (deli .eq. O . and. dely .eq. O) then 
f=0.5 
else if(ynint(sqrt(floaty(deli#delit+tdely#dely))) .eq. 1)then 
f=0.5 
else if(ynint(sqrt(floaty(deli#delitdely#delj))) .eq. 2)then 
#=0. 5 
else 
#=0. 5 
end if 


avdens=0. 0 


do 100 n=0, longest 
scale=0. 5+f 
if (nm .eq. O) then 
avdens=avdens+ro(itn#xinc, ytn*¥yinc, k+ynint(n#zinc) )#Ff 
else if (n .eq. longest) then 
avdens=avdens+ro(i+n#xine, ytn#yinc, kt+ynint(n*#zinc))#0. 5 
else 
avdens=avdenstro(itn#xinc, ytn#yinc, kt+ynint(n¥zinc) ) 
end if 


100 continue 
avdens=avdens/(floaty(longest-1)+scale) ‘!divide by the # of sample pts 
200 continue 


return 
end 
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Appendix 12. Comparison of Heterogeneous Dose Spread Arrays 
Calculated by the Convolution Method and by 
the Monte Carlo Method 
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Table A 
Ai 

=3 ty) =i ) 1 2 3 

) a S 20.12 291 3022 = s 
E = 20102 3081 Obi = = 3 
0 = 1.0g/cm 
hs 3 
1 = .001 .034 .306 .034 -001 ks p = 0.3g/cm 


- -002 -035 281 -035 2-002 - 


2 - 003 . 032 5 Se -032 003 - 
-O001 -004 033 -124 -033 -004 -O01 

Ak 

3 -001 -004 024 -067 ~024 -004 ~O01 
-OO1 -005 023 -061 -023 -004 -001 

4 001 ~005 SOUS) 038 2019 005 O01 
-002 -005 ~OL7 -032 -016 -005 002 

5 002 -005 2013 021 -013 -005 002 
-002 -005 -012 -018 2012 ~005 -002 

6 -001 -004 -008 AONE: -008 -004 O01 


-002 -004 1008" JO -008 -004 -002 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -2.8% 
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Table B 
Ai 

-3 -2 -1 (0) 1 2 3 

) - = 20825. .329 012 = = 
- = 012 ra 27 .012 = e : 
p = 1.0g/cm 
Sele 
1 - .002 .036 2320 .036 .002 = o = 0.2g/cm 


- - 003 -037 - 290 037 - 003 - 


2 -O01 - 003 -038 EOS 038 003 -O01 
-O001 -005 - 036 -135 036 005 -O001 

Ak 

3 -O01 005 028 -081 -028 -005 001 
-002 - 006 -028 -072 -028 - 006 - 002 

4 002 - 006 023 -048 023 -006 002 
-002 - 006 - 020 041 -020 - 006 -002 

5 -002 -005 -016 ~029 -016 -005 002 
-002 -005 -014 -024 -014 -005 002 

6 002 -005 sOLE -018 SOIL 2005 002 


-002 -005 ZOLO -015 -010 -005 002 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -1.3% 
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Table C 
Ai 
=3 =2 =5) 0 1 2 3 
0 = = .012 = = 
a = OZ — — 
" a -001 .025 .239 .025 .001 “= 
= .001 .025 .239 .025 .001 am 
op = 1.0g/cm 
2 = .002 .023 .086 .023 .002 ss p = 0.3g/cm 
001 .003 .021 .074 .021 .003 001 
Ak 
3 = .003 .018 .044 .018 .003 = 
.001 .004 .016 .032 .016 .004 .001 
4 001 .003 .013 .023 .013 .003 001 
.001 .004 .010 .016 .010 .004 .001 
5 001 .003 .007 5012 .007 .003 001 
.001 .003 .006 .008 .006 .003 001 
6 001 .002 .005 .007 .005 .002 .001 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
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Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -2.1% 
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Table D 


Ai 


-OO1 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -0.9% 
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Table E 
Ai 
-3 -2 -1 0 1 2 3 
0 - - ~O12 329 ALO) P - =- 


2 x 2001) 020" 40787 2020), 001 z 
= 001 .020 .073 020 .001 4 
Ak 
3 - FOOT 22008) )020) 2008) 001 2 
2 00% 068  .020 {008 .001 = 
0 = 1.0g/cm> 
Sms eee eo Ee ee ee een ee 
4 - S001. oeawer 20075 «004. .001 5 p = 0.3g/cm” 
“ 001 3003 006 .004 001 = 
5 = 00% 2003-12004" .2003° 1.001 = 
= 0Ol | (3002 "00s. {002° 2001 = 
6 = s00L =50020) 2002 555.002) {00% E: 


- -001 -001 -002 -O01 -O01 - 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -0.2% 
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Table F 
Ai 
=3 -2 -1 0 iL 2 3) 
) 
1 
2 - -001 ~007 -026 ~007 OOM - 
- - - 007 -029 -007 = - 
Ak 
3 - - -003 -008 -003 = - 
- - -003 -009 - 003 - - 
4 - - SOON -002 -OO1 - - 
- - -O01 -002 ~OO1 - - 
5 - - - -OO1 = 2 ea 
~ - - -001 - = = 
6 oy _ = - aS — = 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To +4.3% 
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Table G 
Ai 
BS) -2 =) 0 1 2 3 
fo) = = es = 
1 S O01" wank)? ono sold . .001 = 
a Rant eaekOt4 (otis © "0m ~.001 if, ‘ 
p= 0.3g/cm™ 
3 
2 = Obie | CARLO! Sate 480101 | s0Gd = p = 1.0g/cm 
ss = “Oto):. “4050 5040 °.002 zs 
Ak 
3 a “001.4 12006 Wisi, as006: 15001 zo 
= A 2006 —.019) “2006 4001 a 
4 = “O01. 9002 «604. .002 00 z, 
= * .002 .005 .002 = z 
5 ss = O01) 6.061). 1/4001 S bs 
2 = "Galt *s002 "2001 A a 
G = = = = = re a 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To +7.9% 


a 


Ge, 
= Eeth, 
G8. 
. baK. 
a ia. 
= (th a 
s . 
al 
= 


pl GEi, 
“2 Bii. 
Gio tbc. 
Aten es 
800, EIO, 
Fo aie) $10: 
£60. 


syeq5A Sawwg Srokeusnpegono” sost 
lus09go1s3eh Se yo 


iefunie 


o 1veDom 


Tegint e2 andaok z90ql 
i es el ret st 


#004 of ‘aeaaities re 


># 


449 


Table H 
Ai 
as 3 Bs! 0 1 2 3 
0 = = .012 .012 2 = 
s = .012 .O1l 28 - 
1 = L001y . 025). 2230) 025) 2001 2 
s 00% "02m 289) " 025) 2001 - : 
9 = 1.0g/cm 
2 = 002) .024. = 2083, (0245 .002 = 
E 006  ..008) 0% 023: .007 = 0 = 0.001g/em> 
Ak 
[2 age een ee ee eerie ee eee See 
= 1.0g/cm> 
3 “OOS trrOte 62043 OTs". 008 = p = 1.0g/cm 


-O01 -004 OLT 034 -016 004 O01 


4 - -002 -007 -012 -007 002 - 
- -002 007 -012 ~007 -002 - 


5 - -O01 - 003 -004 - 003 O01 - 
= -O01 003 -004 - 003 O01 - 


6 = - -OO1 001 ~001 - ~ 
~ - -OO1 -OO1 ~OO1 - - 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -1.3% 
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Table I 
Ai 
ZS =o =i ) 1 3 3 
0 es &, FOrz 2S 7OC Tee ONe - zs 
“4 S O12 6329 -Oll = 
1 S .001 .025 $23 000 025 -001 = 
a -001 .025 sos? * 5025 .001 - 
D = .002 .024 .093 .024 .002 2 
.001 .002 .020 .072 .020 .002 = . 
2) 
Ne 0 = 1.0g/cm 
3 .001 .004 .022 .057 .022 .004 .001 
.003 .005 .016 5033 .016 .005 .004 
p = 0.001g/cm” 
4 .001 .004 .014 .027 .014 .004 .001 
.001 .004 .012 ROOM O12 .004 .001 
[an ae Pri cor a ge TR ee Ler ee te ee 
3 
5 = .002 .005 .008 .005 .002 uo o = 1.0g/cm 
.001 .002 .005 .008 .005 .002 .001 
6 = .001 .002 .003 .002 .001 = 
* .001 .002 .003 .002 .001 s 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -1.5% 
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Table A’ 
Ai 
3 ~2 a1 0 1 2 3 
0 = " £012 .012 a * 
: s .012 .O11 2 2 3 
0 = 1.0g/cm 
1 fe 100 77033 4283) 4033 .002 ~ 9 = 0.3g/cm” 


2 - - 003 -03 1 -118 03 1 - 003 - 
-001 -004 - 033 -124 033 -004 -OO1 

Ak 

3 O01 004 023 -059 023 -004 -O01 
O01 -005 023 - O61 - 023 -004 O01 

4 -O0O1 -005 -018 032 -018 -005 -O0O1 
-002 -005 ~O17 -032 -016 ~005 -002 

5 -002 -005 013 -018 -013 -005 -002 
-002 -005 -012 -018 ~012 -005 -002 

6 -O001 -004 -008 sOL -008 004 -O01 


-002 - 004 -008 s@ Lil -008 004 -002 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To 0.6% 
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Za 
Table B 
Ai 

-3 -2 -1 0 1 2 3 

0) - - KORA 0329 AON = = 
- - SOR SIMA ~012 = = 3 
0 = 1.0g/cm 
(oy (O74 /om> 

1 - -002 -035 -290 ZOS0 -002 - “28 


2 -001 003 035 131 -035 - 003 O01 
-001 -005 - 036 also 036 -005 -001 

Ak 

3 -001 -005 -027 -068 -027 ~005 -001 
-002 - 006 -028 -072 -028 - 006 -002 

4 -002 - 006 022 -040 ~022 -006 002 
-002 -006 -020 041 -020 - 006 002 

5 002 -005 -015 -024 -O15 -005 -002 
002 005 -014 -024 -014 -005 -002 

6 -002 -005 -O1l -015 oOL1 -005 002 


-002 -005 OeE -015 -010 -005 -002 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved to 0.3% 
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Table C 
Ai 

-3 -2 oh 0 1 2 3 
e EY ajay 12329) |b 02 2 “ 
a Tai am) ace onl ‘ i 
- O01 ~025 APES) BOs) -O01 = 
SO soo, cous, «eae, ove “2008 - 

0 = 1.0g/cm” 

3 

2X" goo) 2023, «08a 202, ‘Aco? Z o = 0.3g/em 
-OO1 ~003 Ogi 074 O21 ~003 O01 
- -003 018 041 -018 .003 x 
O01 ~004 OLS -032 O16 004 ALOKO)L 
SOO 7003) # Oda aOeth. > 2an2e* 003)" =F 001 
O01 004 -010 016 -010 004 O01 
Wool | .2003) 007, seOogn scom™ = .00s, | 200m 
“00 00s te0e )f0cee S00e. 9003. 2001 
001. .002  .005 .006 .005 .002 .001 
"001 «.002.—S««w004.—Ss«004.—Ss«w004'—s=s«w002~—SC«sw 001 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -1.6% 
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Table D* 
Ai 
-3 -2 at 0 1 2 3 
0 = a .012 |.329 ] .012 - a 


1 iz 001 .025 .239 .025 001 = 
2 .001 .025 .238 .025 -001 = 
2 = .001 .020 .073 .020 .001 = 
001 O01 .020 073 .020 .001 4 
3 
Ak 0 = 1.0g/cm 
ee ——————————— 
3 
3 = .001 O11 .026 SORE 001 = po = 0.3g/cm 
= .002 .010 .022 .010 002 2 
4 = .002 .007 .012 .007 .002 = 


5 ~ O01 -004 - 006 004 ~O01 - 
-O01 -O01 - 003 004 -003 -001 - 
6 -001 -001 003 003 003 -OO1 


-OO1 -001 -002 -002 ~002 -001 -001 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -0.7% 


Te Oe ae Re ee eI H 
ty 


rom. TOs. 905 


iS rt 


eyes PAPICE S50 sucens 
wer eeet ervenopret48 427 20. agi ré tom 


455 


/ 
Table E 
Ai 
-3 -2 -1 ) 1 2 3 
0 - - ~O12 - = 
- - -012 - - 
1 - -OO1 ~025 «299 -025 OO - 
- -OO1 ~025 1239 -025 -O01 - 
2 ~ -OO1 -020 -073 -020 -001 - 
- -O01 -020 OMS - 020 OOM - 
Ak 
3 - -OO1 -008 -020 -008 -O01 - 
- -001 -008 -020 -008 -O001 - 5 
p= 1.0g/cm> 
cence nn ee 
4 - -OO01 -004 -007 -004 -O01 = p= 0.3g/em” 
- -O01 - 003 - 006 -004 OO - 
5 - -OO1 -003 -004 -003 -OO1 - 
- -OO1 -002 - 003 -002 OO - 
6 - -001 - 002 -002 -002 -OO1L - 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -0.2% 
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ae 
Table F 
Ai 

-3 ~2 =a 0 1 2 3 

fe) = = .006 | .143 | .006 " ie 

= = .006 |.145 | .006 . E 
p= 0.3g/cm” 
1 2 y Ola. 098) Oi a = 0 = 1.0g/cm” 


2 - O01 007 028 007 O01 - 
- - 007 -029 ~007 - - 

Ak 

3 - = -004 097 -004 ~ - 
= - - 003 -091 - 003 - - 

4 - - -001 003 -001 - - 
- - -OO1 002 O01 - - 

5 = = = -O01 - - = 
> a = -O01 = ad = 

6 2 & = = 2 = Bs 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -1.3% 
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Table G 
Ai 
=i = -1 0 il 2 3 
0 a 2 .006 . a 
2 - .006 u é 
1 . Toot. sOl4. 21S’ 014) 7).001 e 
3 (Ool .01  otie? ome 001 S 5 
p = 0.3g/an 
3 
2 = rOGL ~~ O10~=S 043 LO. 001 a posed Og/cn 
3 i “O10! 650! .0%0; .001 S 
Ak 
3 = .001 .006 .014 .006 .OO1 2 
s fs “G06, 010? i006 00% 3 
4 a 001. .002. 004 002° .001 z 
5 $ “002 .005 .002 : i 
5 b f 001 .001 .001 2 2 
a ? "001 .002 # .001 M 3 
6 = c = = = 2 = 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To +5.5% 
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A 
Table H 
Ai 
=o =2 -1 0) 1 2) 3 
(0) - - -O12 1329 sO - = 


1 - -001 -025 239 -025 O01 ~ 


S .001 .025 .239 .025 .001 = : 
9 = 1.0g/cm 
2 ES .002 e024 .089 .024 .002 a 3 
a .006 .023 074 .023 .007 E o = 0.001g/cm 
Ak 
ia 3 
3 a .003 ‘O17 $040 “017 *.003 us p = 1.0g/cm 
.001 .004 SOlT .034 .016 .004 .001 
4 4, .002 .006 .O11 .006 .002 - 


- 002 .007 -012 -007 -002 - 


5 - 001 - 003 004 003 -OO1 - 
- O01 - 003 -004 - 003 001 - 


6 - = -001 -001 ~001 - - 
- = O01 O01 O01 - - 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -0.6% 
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ya 
Table I 
Ai 
-3 -2 -1 (0) 1 2 3 
0) - - OZ oa Bs 
- - SOZ = Ss 
1 - -001 O20 0239 2025 -001 = 


2 a .002 .024 .093 .024 002 = 
-001 .002 .020 4072 .020 .002 s : 
Ak 9 = 1.0g/cm 
3 001 .004 3022 .052 A022 .004 -001 
. 003 .005 016 .033 .016 .005 .004 
e = 0.001g/cm> 
4 .001 .004 .014 .024 .014 .004 .001 
.001 .004 2012 .019 BO .004 .001 
3 
5 = 002 005 007 005 -002 2 eo = 1.0g/cm 
001 .002 .005 .008 .005 -002 .001 
6 = .001 .002 .002 .002 .001 = 


- -001 -002 -003 -002 -001 - 


Upper Number Is Calculated From Homogeneous Dose Spread Arrays. 
Lower Number Is From Monte Carlo Simulation Of The Heterogeneous Phantom. 


Energy Conserved To -0.8% 
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PUBLICATIONS: 


The following was a critique solicited by the Cluff Lake Board of 
Inquiry chaired by Justice Bayda and delivered to the Inquiry in 
Regina, Saskatchewan, July 1977: 


1. Review of "Appendix C, Radiological Health and Safety", 
Caplan, H., Mackie, T.R. 


The following were published under the auspices of the Atomic Energy 
Control Board of Canada under the aegis of the Federal-Provincial 
Task Force on Radioactivity: 


1. Mackie, T.R., ''Reducing Airborne Radiation in Dwellings 
using Make-Up Air Ventilation'', March 1979. 


2 Curly, Re, Mackie; TR, Haubrich, E.J.J., "A Study of Air 
Exchange Rates in Dwellings using SF, as a Tracer Gas 
May 1980. 


SueGrillys Ro, Mackie, Wena: yar Make-Up Units for Reducing 
Radiological Levels'', May 1980. 


The following article is reprinted in Appendix 9: 


1. Mackie, T.R., Scrimger, J.W., "Contamination of a 15 MV 
Photon Beam by Electrons and Scattered Photons'', Radiology, 
144, 403-409, July 1982. 


The following articles appeared in the proceedings of ''The Eighth 
International Conference on the Use of Computers in Radiation Therapy", 
July 1984: 


1. Mackie, T.R., Scrimger, J.W., "Computing Radiation Dose for 
High Energy X-rays using a Convolution Method". 


Battista, J.J Mackie, lho: El-Khatib, E.; Serimger, J.W., 
"Tung Dose Corrections for 6 MV and 15 MV X-rays: Anomolies". 


Be Mackie, [eho pabtistase)<) «.00 4 Macroscopic Monte Carlo 
Method for Electron Beam Dose Calculations: A Proposal". 


The following articles have been submitted to the journal, "Medical 
PHySics::: 


1. Mackie, T.R., Scrimger, J.W., Battista, Jada a COMVOLUCTON 
Method of Calculating Dose for 15 MV X-rays". 


2. Mackie, T.R., El-Khatib, E., Battista, Users SOTA Ve ew elicg, 
Van Dyk, J., Cunningham, J.R., ''Lung Dose Corrections for 
6 MV and 15 MV X-rays". 
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